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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

PERIODIC  OPTIMIZATION  OF  BIOLOGICAL  SYSTEMS 

by 

Eva-Maria  Abulesz 
August  1 988 

Chairman:  Dr.  Hong  H.  Lee 

Major  Department:  Chemical  Engineering 

In  bioprocessing  and  medicine  periodic  conditions  frequently  prove 
to  be  more  successful  in  achieving  a certain  objective  than  steady-state 
conditions.  Often  this  is  due  to  the  fact  that  cells  need  some  time  to 
adapt  to  new  conditions,  when  their  environment  changes.  If  one  can  use 
this  lag-time  to  aid  in  reaching  the  goal  of  interest,  periodic 
operating  conditions  are  highly  desirable. 

In  this  dissertation  methods  of  optimal  control  theory  are  applied 
to  establish  the  optimal  periodic  conditions  in  two  major  areas: 
fermentation  technology  and  cancer  chemotherapy. 

The  method  of  ir-criterion  is  used  to  determine  whether  periodic 
variation  of  the  dilution  rate  can  enhance  the  performance  of  continuous 
fermentation  processes.  It  is  found  that  the  presence  of  a time-delay 
in  the  dynamic  response  of  the  chemostat  renders  a periodic  operation  of 
bioreactors  which  results  in  enhanced  biomass  producti vity . Also 
employing  Williams’s  structured  model  it  is  shown  that  cycling  can 
improve  the  average  protein  productivity. 
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Experiments  with  S.  oerevisiae  are  described  which  aim  at  verifying 
these  predictions.  The  dynamic  model  found  to  describe  step  changes  in 
dilution  rate  incorporates  a lag-time  during  which  the  culture  adapts  to 
the  new  environmental  conditions.  It  was  confirmed  that  for  certain 
cycling  frequencies  periodic  reactor  operation  results  in  enhanced 
biomass  productivity. 

The  ir-criterion  is  then  applied  to  cancer  chemotherapy.  A 
procedure  is  proposed  which  allows  to  determine  the  optimal  treatment 
regimen  given  a certain  patient  and  drug.  As  a result  it  can  be  seen 

whether  the  drug  under  consideration  will  have  the  desired  effect  as 
well  as  which  form  of  treatment  would  be  preferable:  cycling  (on/ off) 

or  continuous.  In  the  event  cyclic  treatment  is  superior,  the  optimal 
treatment  regimen  is  established. 

To  assess  the  effect  a series  of  injections  of  drug  has  on  the 
normal  and  malignant  cell  population,  and  on  the  performance  measure,  a 
new  method,  "periodic  impulse  forcing  of  nonlinear  systems,"  is 
developed.  The  method  is  based  on  the  Carleman  linearization  procedure 
and  allows  the  determination  of  the  state  equations  and  the  performance 
measure  as  a function  of  the  forcing  period.  It  is  subsequently  applied 
to  a model  for  cancer  chemotherapy. 
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CHAPTER  I 

GENERAL  INTRODUCTION 
Motivation  and  Background 

In  the  area  of  biochemical  engineering,  although  substantial  re- 
search effort  has  been  invested  in  the  optimization  of  batch  and  semi- 
batch processes  [1-6],  very  little  has  been  done  in  the  optimization  of 
continuous  bioreactors  [7].  Over  the  past  few  years  it  has  been  real- 
ized that  steady-state  operation  of  chemical  reactors  is  not  necessarily 
the  optimum  type  of  operation.  If  some  manipulated  process  variables 
are  allowed  to  vary  with  time,  very  often  superior  productivity,  para- 
metric insensitivity  or  even  a change  of  the  culture  fate  is  possible 
when  compared  to  steady-state  operation.  A significant  amount  of  work 
has  been  undertaken  [8-14]  to  establish  the  periodic  operations  that 
maximize  the  yields  and/or  selecti vities  of  chemical  reactors.  Some 
workers  [15-20]  have  tried  to  assess  the  effect  of  cycling  behavior  of  a 
chemostat  culture  experimentally,  but  it  is  hard  to  draw  any  general 
conclusions  3ince  frequently  undefined  growth  media  (e.g.,  molasses)  or 
insufficiently  defined  quantities  were  used. 

The  quality  of  any  optimization  depends  to  a great  extent  on  how 
well  the  mathematical  model  at  hand  describes  the  process  in  question. 
In  general  two  types  of  models  are  used  to  describe  microbial  growth 
processes:  "unstructured  models"  and  "structured  models."  Unstructured 

models  have  inherent  the  assumption  of  balanced  growth.  This  indicates 
that  every  single  part  of  the  cell  grows  at  the  same  rate  during  a 
period  of  balanced  growth.  This  situation  may  be  viewed  as  the  state 
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reached  after  all  the  transients  have  died  out  and  the  cell  has  adapted 
to  a constant  environment,  essentially  at  steady-state.  For  this  very 
reason  unstructured  models  succeed  in  describing  steady-state  situa- 
tions, but  they  are  in  general  insufficient  in  accurately  portraying 
dynamic  systems.  For  such  systems  the  assumption  of  balanced  growth  has 
to  be  relaxed.  This  can  be  done  in  two  ways.  The  first  possibility  is 
to  include  a time-delay  in  the  mathematical  description  of  the  system. 
Introducing  a time-delay  increases  the  dimensionality  of  the  mathemati- 
cal model  but  it  allows  a better  description  of  dynamic  situations. 
Moreover  it  is  also  physiologically  meaningful  because  it  essentially 
assumes  that  the  cell  does  not  react  to  a change  in  environmental 
conditions  by  an  instantaneous  adaptation  of  its  growth  rate. 

The  alternate  way  to  relax  the  assumption  of  balanced  growth  is  the 
formulation  of  a structured  model.  Such  models  break  down  the  cell  in 
separate  compartments,  as  many  as  needed,  but  at  least  two,  and  allow 
for  each  compartment  to  grow  at  its  own  rate.  Structured  models  are 
superior  to  unstructured  ones  in  describing  dynamic  situations.  It  is 
not  very  meaningful , though,  to  use  models  of  too  high  dimensionality 
because  aside  from  their  mathematical  complexity,  many  of  the  state 
variables  in  such  models  (nucleic  acids,  metabolites,  growth  factors) 
are  not  readily  measurable,  especially  on  line,  something  that  makes  the 
determination  of  system  parameters  practically  impossible. 

Based  on  these  various  types  of  mathematical  models  methods,  of 
optimal  eontrol  theory  will  be  used  to  detect  the  desirability  of  peri- 
odic reactor  operation.  Moreover , if  periodic  operation  is  the  opera- 
tion of  choice,  an  estimate  of  the  optimal  cycling  frequency  will  be 
obtained.  Experiments  will  be  conducted  to  assess  the  effect  of  cyclic 
reactor  operation  on  yeast  in  a chemostat  culture. 
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Just  as  the  growth  of  micro-organisms  can  be  stimulated  by  a cer- 
tain "feeding  schedule,"  the  growth  of  tumor  cells  can  be  minimized  by  a 
certain  "treatment  schedule."  Unfortunately  the  situation  here  is  much 
more  complex  than  in  the  case  of  uni-cellular  organisms.  So  far  the 
success  of  chemotherapy  and  radiation  therapy  has  been  limited.  Surgery 
proves  ineffective  in  about  half  of  all  cases  because  the  tumor  was 
discovered  too  late  and  metastases  have  already  spread  throughout  the 
body.  So  it  is  not  surprising  that  up  to  400,000  people  die  from  cancer 
every  year  in  the  U.S.  alone. 

Healthy  mammalian  cell3  are  under  strict  homeostatic  control.  New 
cell  formation  is  turned  on  or  off  according  to  the  body's  need.  In 
cancer  cells,  however,  the  control  mechanism  breaks  down  and  uncon- 
trolled growth  occurs.  Therefore  a cancer  cell  has  a selective 
advantage  over  a benign,  controlled  cell.  Cancer  cells  in  general  do 
not  grow  faster  than  normal  body  cells.  Cancer  cells  found  even  in 
constantly  proliferating  tissues  like  epitheliiim  or  bone  marrow  have 
longer  cell-cycle  times  than  their  benign  stem  cell.  What  makes  cancer 
such  a dreadful  disease  and  its  treatment  extremely  difficult  is  its 
invasiveness  of  other  body  compartments,  making  localized  surgery  inef- 
fective. Recent  alternatives  like  monoclonal  antibodies  are  still  in 
their  infant  stages  and  so  chemotherapy,  though  not  an  optimal  form  of 
treatment,  has  to  be  improved  in  any  way  possible.  The  approach 
followed  in  this  dissertation  is  based  on  the  mathematical  modeling  of 
normal  and  cancer  cells.  The  kinetic  differences  between  those  cell 
types  are  used  to  optimize  the  treatment  schedule.  Recall  that 
unstructured  models  in  general  fail  to  predict  the  dynamic  of  even  a 
unicellular  organism.  That  is  even  more  so  for  such  a complex  structure 
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as  the  human  body.  Because  it  is  not  an  ideal  CSTR  (continuous  stirred 
tank  reactor)  with  complete  mixing  it  is  frequently  necessary  to 
introduce  models  comparable  to  the  structured  models  introduced  above. 
In  pharmacokinetic  modeling  the  body  is,  if  needed,  divided  into  various 
compartments  each  exhibiting  similar  kinetic  properties.  These 
compartments  may  or  may  not  comprise  physical  entities,  such  as 
organs.  But,  of  course,  the  model  should  be  kept  as  simple  as  possible, 
describing  all  the  known  facts  and  able  to  support  the  intended 
purpose.  For  this  reason  this  work  is  limited  to  one-compartment 
pharmacokinetic  models  and  cycle  non-specific  drugs  as  a first  stepping 
stone  towards  predictable  chemotherapy. 

Dissertation  Outline 

In  this  dissertation , biological  processes  like  continuous 
f ermentations  and  drug  therapy  will  be  described  by  a mathematical 
model.  Using  two  methods  of  periodic  optimal  control,  the  optimal 
periodic  input  (or  the  optimal  steady  state  operation,  if  it  is  found  to 
be  superior)  is  determined.  Experiments  are  described  which  assess  how 
periodic  input  actually  influences  the  chemostat  culture. 

In  Chapter  II  the  optimization  theory  used  in  this  thesis  will  be 
presented.  Following  an  overview  of  relevant  optimal  control  theory  and 
the  Tr-criter ion  (II.  1.)  a new  method  will  be  introduced  (II. 2.), 
concerned  with  periodic  impulse  forcing  of  nonlinear  systems.  It  allows 
the  explicit  evaluation  of  the  performance  measure  of  such  a system  as  a 
function  of  the  cycling  frequency. 

Chapter  III  deals  with  the  application  of  periodic  optimal  control 
to  continuous  fermentations.  In  the  first  part  (III.1)  the  desirability 
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of  periodic  operation  is  evaluated  for  unstructured  models,  unstructured 
models  with  time-delay  and  structured  models.  In  the  second  part 
experimental  results  on  the  behavior  of  Saccharomyces 
cer e vi si ae  (yeast)  in  a chemostat  culture  under  periodic  square  wave 
forcing  are  presented. 

In  Chapter  IV,  optimal  control  theory  is  applied  to  cancer  chemo- 
therapy. In  the  first  part  of  this  chapter  (IV. 1.)  it  is  shown  how  one 
can  determine  whether  periodic  or  continuous  treatment  of  cancer  renders 
better  results  and  how  appropriate  the  drug(s)  used  is  (are)  in  a parti- 
cular situation. 

The  new  method  described  in  Section  II. 2.  will  be  applied  to  drug 
therapy  in  Section  (IV. 2.)  to  assess  the  effect  of  repeated  injections 
(periodic  impulse  forcing). 

Chapter  V finally  draws  general  conclusions  obtained  from  this 
study  of  periodic  optimization  of  biological  systems. 


CHAPTER  II 

METHODS  OF  OPTIMAL  CONTROL  THEORY 
AND  PERIODIC  OPTIMIZATION 

II.  1.  Optimal  Control  Theory  and 
the  n-Criterion:  An  Overview 


Introduction 

To  apply  optimal  control  theory  to  an  optimization  problem  requires 
the  following: 

• a mathematical  model  of  the  process  to  be  controlled 

• the  knowledge  of  all  constraints 

• a known  objective,  the  perf ormance  criterion 

The  classical  control  problem  is  the  one  where  an  object  shall  be 
moved  from  a state  A to  a state  B.  The  object  might  be  a car,  a rocket 
or  the  vector  oi  state  variables  for  a chemical  reaction;  the  states 
might  be  points  in  space  or  the  steady  states  of  a reaction  system. 

The  state  of  the  car,  for  example,  can  be  described  by  its  position 
and  i s^s  velocity  (state  variables).  It  can  be  controlled  by  varying  the 
acceleration  and  the  brake  deceleration  (control  variables)  of  the 
vehicle.  The  task  to  perform  can  be  to  move  the  car  from  point  A to 
point  B on  a straight  line.  Its  motion  can  be  described  by  the 
equations  of  classical  mechanics  (mathematical  model  of  the  process). 

In  general  three  kinds  of  constraints  can  be  distinguished:  state 

constraints,  which  define  the  values  of  the  state  variables  at  points  in 
space  and  time — often  called  "boundary  conditions"  (e.g.  at  the  initial 
and  the  final  time  the  car  is  stationary);  control  constraints,  which 
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set  limits  on  the  admissible  values  for  the  control  variables  (e.g.  the 
highest  possible  acceleration  is  ...);  and  additional  constraints  (e.g. 
the  amount  of  gas  in  the  tank). 

Knowing  the  model,  the  constraints  and  the  task  to  be  performed, 
the  question  arises  what  is  the  "best,"  the  "optimal"  way  to  perform 
this  task?  What  are  the  criteria,  that  make  one  trajectory  (history  of 
state  variables)  better  than  the  other?  If  the  objective  is  to  move  the 
car  from  A to  B in  the  least  amount  of  time,  then  the  fastest  admissible 
trajectory  is  better  than  any  slower  one.  If  the  objective  is  to  use 
the  least  amount  of  gas  on  the  way,  then  a slower  trajectory  might  be 
optimal.  The  objective  is  reflected  in  the  "performance  measure"  or  the 
"performance  criterion." 

In  the  same  way  as  this  problem  of  classical  mechanics  can  be 
solved,  biological  systems  can  be  optimized.  Whether  the  objective  is 
the  fastest  growth  of  a microorganism , the  highest  productivity  of  an 
antibiotic  or  the  maximum  kill  of  tumor  cells,  the  principles  stated 
above  are  still  valid. 

In  the  next  section  the  optimization  problem  will  be  formulated 
mathematically . 

Mathematical  Formulation  of  the  Optimization  Problem 

Let  the  process  to  be  optimized  be  defined  by  a set  of  ordinary 
differential  equations  (although  other  descriptions  such  as  partial 
erential  equations  are  also  allowable),  which  adequately  describes 
the  response  of  the  system  to  all  possible  inputs.  Assuming  that  f has 
continuous  partial  derivatives  in  x and  u, 
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x(t)  = f (x(t),  u(t),  t) 


(2-1) 


where 


x(t) 


x1  (t) 
x2(t) 


xr(t) 


(2-2) 


is  the  vector  of  state  variables  or  state  vector  and 


u(t)  = 


u1  (t) 
u1  (t) 


us(t) 


(2-3) 


is  the  vector  of  control  variables  or  control  vector. 

In  most  cases  the  objective  of  the  process  can  be  fulfilled  by  a 
multitude  of  different  control  vectors,  called  the  admissible  con- 
trols. The  purpose  of  applying  optimal  control  theory  is  to  determine 
the  admissible  control  which  maximizes  or  minimizes  the  performance 
criterion  J which  is  given  by 

V 

J = / <p  ( x ( t ) , u(t),  t )dt  + ip  (x(t  ) , t_)  (2-4 ) 

4-  “it 

0 

where  tQ  and  t^-  are  the  initial  and  the  final  time,  respectively,  and 
4>  and  \p  are  scalar  functions,  possessing  continuous  partial  derivations 
in  x and  u.  The  first  term  in  equation  (2-4)  depends  on  the  state 

trajectory,  which  in  turn  is  determined  by  the  control  history  and  the 
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initial  conditions,  over  the  time  interval  At  = tf  - t , whereas  the 
second  term  is  a function  of  the  final  state  only. 


The  Pontryagin  Maximum  Principle 

An  optimal  control  u^j.  maximizes  the  performance  criterion  J if 


AJ  = J(uQpt)  “ J(u)  - 0 (2-5) 

Any  control  u.  (if  sufficiently  close  to  u^^)  can  be  expressed  as 


u 


u 

-opt 


+ 5u 


(2-6) 


so  that  the  necessary  condition  to  minimize  the  functional  J becomes 


5 J (u  . , 6u)  > 0 
-opt 


(2-7a) 


if  Jiopt  lles  on  the  boundary  of  admissible  values  for  the  control  u at 
any  time  during  the  interval  At,  or 

'5J(-opt’  5-)  = 0 (2-7b) 

-liopt  lies  withiR  the  boundaries  in  the  time  interval  At. 

Recall  now  the  functional  to  minimize  as  given  by  equation  (2-4). 
Rewriting  the  performance  measure  as 

t„ 

t 

J = J U(x(t),  u(t),  t)  + lT(t)f[x(t),  u ( t ) , t)  - X ( t ) ] 

° (2-8) 

- i|»(x(tf),  tf) 
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where  X is  the  vector  of  adjoint  variables  or  Lagrange  Multipliers,  and 
defining  the  Hamiltonian 


H(x(t ) , u(t),  X (t ) , t)  = <j>  (x  (t ) , u(t),  t) 
+ XT(t)f (x(t),  u(t),  t) 


(2-9) 


it  is  obvious  that 


J = I [ H(x  (t ) , u(t),  t)  - XT(t)x(t)]dt  + ip  ( x ( t f ) , t ) (2-10) 

t 


which  can  be  integrated  by  part  to  yield 


J = / [H(x(t),  u(t),  t)  + XTx(t)]dt  + \p  (x  (t  „ ) , t 

t " " " f f 

o 


) - XTx 


f (2-11) 
o 


It  can  be  shown  [21 ] that  to  minimize  J it  is  necessary  that 


H<5opt(t>’  Sopt(t)  * iopt(t)'  U S 

(2-12) 

H<2opt(t)’  Sopt(t)'  V(Ul  ° 

Equation  (2-12)  is  the  essence  of  the  Por.tryagin  Minimum  Principle.  It 
states  that  to  maximize  the  performance  criterion  J an  optimal  control 
Uopt  has  to  maximize  the  Hamiltonian  H. 
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Summarizing,  for  a process  described  by  equations  (2-1)  - (2-3)  and 
a performance  measure  as  given  by  equation  (2-4)  to  be  optimized  by  an 
optimal  control  vector  u^^ , the  Hamiltonian  is  defined  as  in  equation 
(2-9).  The  following  conditions  must  hold  for  t e At: 


(2-13) 


9H 

3A 


(2-14) 


A 


3H 

3x 


(2-15) 


Equations  (2-13)  - (2-15)  define  a system  of  equations  in  x and  A.  A 
two-point  boundary  problem  has  to  be  solved.  The  boundary  conditions 
are  provided  by  the  transversality  conditions 

" i)  = 0 (2-16) 

for  t = tQ  and  tf. 


The  ir-criterion 

As  discussed  in  Chapter  1 in  many  cases  periodic  operation  proves 
to  be  very  desirable.  The  problem  of  determining  the  optimum  periodic 
operation  was  addressed  by  Bittanti  et  al . [22]  who  developed  the  so 
called  ir-criterion  for  optimality.  This  criterion  states  a sufficient 
condition  for  improvement  of  the  performance  of  a periodic  system  via 
cycling.  The  method  In  addition  provides  an  estimate  of  the  optimum 
cycling  frequency  for  the  particular  system  of  interest. 
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This  section  sets  the  mathematical  framework  for  studying  periodic 
optimization  problems  by  applying  the  u-cr iter ion. 

Recall  the  model  equations  as  given  by  equation  (2-1)  and  let 

l (t)  = h(x(t) , t)  (2-17) 

be  the  equation  relating  the  measured  variables  y to  the  state  variables 
x.  Assume  that  there  exists  a steady-state  solution  corresponding  to 
u = u°  for  which  £ = _x°  and  y_  = y°,  at  which  a performance  measure 
J = <(>(u,  y)  is  maximized.  We  seek  a periodic  operation  of  the  form 

u(t)  = u ( t + t)  (2-13) 

that  would  lead  to  an  average  performance  measure 

J(u(*);  T)  = “ Jq  y)dt  (2-19) 

which  is  higher  than  J°  = JQpfc  = <J>(uQpt , y^). 

In  most  cases  it  is  a good  assumption  [23]  that  a control  policy  of 
the  form  (2-18)  leads  to  a periodic  solution  of  the  same  period;  i.e. 

x(t  + t)  = x(t)  (2-20) 

The  optimization  is  in  general  subject  to  equality  constraints  of  the 
form 


7 Jo  ’i)dt  = 0 


(2-21) 
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and  inequality  constraints  of  the  form 

1 r x 

7 J0  w(Z»u)dt  * 0 (2-22) 

Assuming  that  f(x,u),  h(x),  <t>(jr,u),  y(y,u)  and  w(y  ,u ) are  twice 
differentiable  at  the  optimum  steady-state,  one  can  define  the  Hamilton- 
ian 

H(x,u,A,jj)  = <)>(h(x),  u)  + ATf(x,  u)  + yTv(h(x),  u)  (2-23) 

X and  y being  Lagrange  multipliers. 

Let 

ir(o))  = G (-jai)PG(ja))  + Q G(jio)  + G (-jto)Q  + R (2-24) 

where 


G(p)  = (pi  - A)  ^ (2-25) 
and  A = fx(x°,  u°)  (2-26) 

B = S°)  (2-27) 
P = Hxx(-°’  -°’  -0’  (2-28) 
Q = Hxu(-°’  -°>  ^°>  (2-29) 
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R = Huu(-°’  -°*  -°'  y0)  (2-30) 

Then  the  TT-criterion  states  that  if  the  (n  x n)  Hermitian  matrix  ir(o>)  is 
partially  positive  for  some  values  of  u > 0,  then  there  exist  to  > 0 such 
that  ir(u))  is  not  negative  definite  [22].  For  single  input-single  output 
systems  of  the  type  concerned  with  this  work,  ir(u)  is  simply  a scalar 
function  of  the  frequency  w.  For  the  optimization  problems  of  interest 
it  then  suffices  to  have  tt  > 0 for  some  cycling  frequencies  u.  This 
will  guarantee  that  operating  at  such  cycling  frequencies  will  result  in 
superior  bioreactor  performance. 

II.  2.  Periodic  Impulse  Forcing  of  Nonlinear  Systems-- 

A New  Method  ' 


Introduction 

The  ir-criterion,  as  introduced  in  Section  (II. 1.),  gives  only 
sufficient  conditions  for  the  improvement  of  the  performance  of  a 
certain  system  by  periodic  operation.  Moreover,  it  is  strictly  valid 
only  locally  and  does  not  give  any  information  about  the  optimal  wave 
form  or  amplitude.  To  obtain  this  information,  tedious  and  time 

consuming  integrations  have  to  be  performed.  A new  procedure  is  needed 
to  obtain  a direct  relationship  between  the  performance  criterion,  the 
forcing  frequency  and  the  system  parameters. 

In  the  practice  of  periodic  operation  the  two  most  easily  implemen- 
table  wave  forms  are  impulse-  and  pulse-forcing.  For  the  case  of  pulse 
forcing,  Lyberatos  and  Svoronos  [24]  utilized  the  Carleman  linearization 
procedure  to  develop  a method  which  allows  the  determination  of  the 
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performance  measure  as  a function  of  the  three  pulse  forcing  parameters, 
which  suffice  to  describe  the  system.  Their  result  is  given  in 
Appendix  A. 

Impulse-forcing  has  widespread  applications  in  many  areas  of  indus- 
try, biology  and  medicine.  Semi  batch-operation  of  a reactor,  periodic 
harvesting  of  a fishculture  in  a pond,  repeated  injections  of  antibio- 
tics into  a ferraentor  or  periodic  injections  of  drugs  in  medicine  are 
just  a few  examples.  Since  this  method  is  easily  implementable , often 
more  economical  and  in  many  cases  improves  the  performance  of  the  sys- 
tem, it  frequently  proves  to  be  the  operation  of  choice. 

The  mathematical  characteri zation  of  this  type  of  operation,  on  the 
other  hand,  often  proves  to  be  rather  difficult  or  impossible.  If  the 
mathematical  model  employed  is  linear,  the  evaluation  of  the  time- 
dependence  of  the  state  variables  is  straightforward.  For  nonlinear 
models,  so  far  only  numerical  methods  have  generally  been  used. 

In  the  following,  the  Carleman  linearization  procedure  is  applied 
to  evaluate  systems  which  undergo  periodic  impulse  forcing.  This  new 
method  is  applicable  to  nonlinear  lumped  parameter  systems  and  perfor- 
mance measures  and  it  allows  the  explicit  evaluation  of  the  time-depen- 
dence of  the  system  as  well  as  the  performance  measure. 

The  next  section  states  the  mathematical  problem,  followed  by  a 
review  of  the  Carleman  linearization  procedure.  Then  the  new  method  is 
introduced  and  illustrated  by  an  example. 
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Statement  of  the  Problem 


Consider  the  nonlinear  system 


l = f(z) 


(2-31) 


with 


f (o)  = o 


(2-32) 


Suppose  a series  of  impulses  is  introduced  with  a period  t which  results 
in  an  instantaneous  addition  of  a disturbance  m to  the  state  vector  y. 
Assuming  that  the  response  is  also  T-periodic,  after  all  transients  have 
died  out,  the  following  relationship  must  hold, 


where  and  y^  denote  the  values  of  the  state  variables  at  the  begin- 
ning and  end  of  the  period  respectively.  The  objective  of  this  investi- 
gation is  to  find  an  analytical  solution  for  the  nonlinear  system  as 
described  by  equations  (2-31)  and  (2-32)  which  makes  it  possible  to 
calculate  explicitly  an  average  performance  measure  J which  can  be 
written  as 


(2-33) 


1 rT 

J = - J P(y)  dt 


(2-34) 


o 


as  a function  of  the  forcing  frequency. 
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The  Carleman  Linearization  Procedure 

The  Carleman  linearization  procedure  is  an  approximation  of  a 
nonlinear  system  by  a linear  system  of  higher  order  [25].  Its  previous 
application  to  nonlinear  system  dynamics  and  control  can  be  found  in 
references  26-35.  The  procedure  will  be  outlined  in  the  following. 

The  function  t_  of  system  (2-31),  assuming  it  is  differentiable  up 
to  order  l at  y = o,  is  expanded  into  a Taylor-series  and  the  monomials 
of  order  up  to  2,  are  introduced  as  new  variables.  When  these  monomials 
are  differentiated  and  only  terms  up  to  order  l are  retained,  a higher 
order  linear  system  in  the  new  variables  is  obtained.  For  example  the 
system 


x 


= (sinx.  )x 


2 


(2- 3 5a) 


(2-35b) 


is  Taylor- expanded  around  zero  into 


( 2-  3 6a ) 


( 2- 3 6b ) 


This  system  has  the  second  order  Carleman  linearization: 
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d_ 

dt 


X1X2 


0 

1 

0 


0 0 
0 0 
0 0 


0 10  0 
0 10  0 

0 0 0 0 

0 10  0 

0 0 10 

0 0 0 2 


(2-37) 


In  general  the  Taylor  series  approximation  of  system  (2-31)  around  y = o 
can  be  written  as 


1 = ^ 
j=1 


A1.J  * 


[j] 


(2-38) 


where  y^  ^ = y_xy_x  ...  xy_(x  denotes  Kronecker  multiplication). 

The  2,th  order  Carleman  linearization  is  the  approximation  of  the 
original  system  by  the  following  linear  system: 


• 

l 

A1,1  A 1 * 2 ‘ Ai,£ 

i[2] 

° A2,1  ‘ A2,£-1 

i H 

• 

• • 

• • 

l:  ' • 1 

1 1.1  / 

(2-39) 

w 


where 


n 


A. 
i .J 


+ 


( 2-40) 
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Solving  Nonlinear  Systems  Undergoing  Impulse  Forcing  by  Using  the  Carle- 
man  Linear 1 zatlon  Procedure  

As  outlined  in  the  previous  section,  the  nonlinear  system  described 
by  equation  (2-31)  can  be  approximated  through  the  Carleman  lineariza- 
tion procedure  to  yield  the  following  linear  system 


! ■ C,» 


(2-41 ) 


where  is  the  i,th  order  Carleman  matrix  and  is  given  by 

• • • • A. 


A1.1  A 1 , 2 

° 


1,2. 


2,1 


(2-42) 


where  is  given  by  equation  (2-40).  This  linear  system  has  the 

solution 


V 

Sf  " 6 *o 


(2-43) 


where  wf  and  are  the  vectors  of  Carleman  variables  at  the  final  and 
initial  time,  respectively,  and  t denotes  the  period.  In  any 
T-periodic  system  the  value  of  w is  given  by 


w = ( (w  - m)  (w  - m)[2]  ...  (w  - m)[£])T 

I “O  - ~o  - -o  - 


(2-44) 


where  _m  is  the  impulse  vector. 
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In  a similar  manner  the  function  P in  the  perf ormance-measure 
(equation  (2-310)  can  be  expanded  into  a Taylor  series  around  0 so  that 
the  performance-measure  takes  the  following  form  in  Carleman  coordi- 
nates: 


1 f ' T - 1 ^ A T 

J = - j r w(t)  dt  = — C [e  1 - I]  w 

^ T 36  -O 


(2-45) 


Theorem.  The  vector  of  state  variables  at  the  beginning  of  a 
period,  Wq  is  given  by 


where 


C,T  4-1  -1 

wQ  = [e  - I A.]  m 
i=0 


(2-46) 


* 

m = 


(-m) 

(-m) 


C2] 


(-m) 


[4] 


(2-47) 


A = I 
o 


(2-48) 


and  for  i > 1 
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-n  x n 


0 . 

- 1 

n x n 


n x n 


E*I  x (-m)Ci] 


A.  = 
1 


-°  £ 

n x n 


(2-49) 


E*I[£~i]x  (-m)Ci]  0 o 

V x n*-1*1  X 


A general  example  is  given  in  Appendix  B.  Here  the  symbol  E*  denotes 
"sum  of  all  possible  orderings;"  for  example 


E [a^b]  = a^b  + a^ba  + aba^  + ba^ 


(2-50) 


Proof.  Substituting  equation  (2-33)  into  equation  (2-43)  gives 


1 

(Z0  - 2) 

* 

io 

% ■ 2?)C2] 

V 

y C2] 
io 

• 

= e 

% - «)ctl 

*o[*] 

Til 

where  (^  - m)L  = (^  - m)  x (y^  - m)  x ...  x (y^  - m) 


Expanding  the  left-hand-side  of  equation  (2-51) 
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+ (-m) 


y0[2]+  E*[lx(-m)]  yQ 


+ (-m) 


[2] 


2^t;l  » j0Ct'1J  * ...  ♦ l "Zlxi-m)'-'  * (-m) 


o 


V 

= e w 


(2-52) 


If  one  defines  the  matrices  as  in  equations  (2-48)  and  (2-49)  and 
substitutes  into  (2-52),  one  obtains 


*0 

,C2] 
j *° 

*o 

*0 

y[2] 

l0 

1*0  1 
42] 

1 H 

+ ai 

> 

io 

+ A3 

: 

ju] 

zo 

■ilz] 

1 

(-m) 

(-s)C21 

CPT 

1 

ZC2] 

*0 

+ + h-2 

l;"J 

+ 

1 

I--)"1 

X, 

= e 

-o 

1 

(2-53) 


Recognizing  the  vector  as  a common  multiplier  and  using  equation 
(2-47),  equation  (2-53)  becomes 


23 


Aw  + 
o-o 


A,w 

1-0 


Vo 


A_w 

3-o 


AH-1-o  +m 


(2-54) 


which  can  be  solved  for  wq  to  yield  equation  (2-46). 

Remark  1 . For  most  systems  it  is  a good  assumption  [22]  that 
inflicting  disturbances  periodically  lead  to  a periodic  solution  of  the 
same  period;  i.e. 


y(t  + x)  = y(t)  (2-55) 

But  as  the  following  example  shows,  there  can  be  systems  or  ranges  of 
admissible  periods  for  which  this  assumption  is  not  valid.  Possible 
outcomes  include  multiperiodicity  or  chaotic  oscillations  [36],  For 
these  systems  the  new  procedure  fails  to  predict  the  occurrance  of 
multiperiodic  or  random  oscillations. 

Remark  2.  For  small  values  of  the  period  x the  matrix  to  premulti- 
ply  _m  in  equation  (2-46)  can  become  nearly  singular.  Because  this 
matrix  needs  to  be  inverted,  numerical  problems  might  result.  In  that 
case  maximum  pivot-strategy  (see  Appendix  C)  should  be  employed. 

Example 

Consider  the  population  growth-model  due  to  Verhulst  [37] 

x = x - x2  (2-56) 

and  assume  that  harvest  takes  place  at  a constant  rate  of  r = — = 3/1 6. 
The  average  performance-measure  is  given  by 

J = — f x(t)  dt 

T 

O 


(2-57) 
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and  the  stable  steady  state  is 


x = 3/4 


( 2-  5 8a ) 


at  which 


J = 3/4 


(2-58b) 


The  steady-state  operation  shall  be  compared  to  periodic  operation  where 
at  each  time-interval  t a part  of  the  population  equal  to  3/16  x will  be 

harvested.  Clearly  in  this  case  m = - -^  x. 

1 6 

A major  concern  of  the  method  presented  is  the  order  of  the  Carle- 
man  approximation  that  needs  to  be  used.  This  is  dependent  on  the 
magnitude  of  the  disturbance,  which  is  proportional  to  x,  the  harvesting 
period.  The  longer  the  period,  the  larger  is  the  disturbance,  the 
higher  is  the  order  of  approximation  needed. 

The  system  described  by  equation  (2-56)  can  be  solved  analytically 
for  x to  yield 


x 


-t 

x e 
o 


1 + x (1-e_t) 
o 


(2-59) 


As  mentioned  in  the  previous  chapter,  it  is  of  importance  to  know, 

whether  the  system  actually  is  x-periodic.  A quick  method  to  check  is 

the  procedure  of  "Poincare  mapping,"  where  x_,  _.1  is  plotted  vs . x 

y i o 9 n 

Figure  2-1  shows  that  the  Verhulst- model  (a=1)  has  three  regions  of  x 
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which  exhibit  different  behavior.  For  m = - 3/1 6t  x-periodicity  is 
predicted  for  o<x<-4.1  and  t > 6.  For  4.1  < t < 6 random  oscil- 
lations take  place.  This  can  be  verified  by  letting  x = x . and 

O t ri  I 

t = t in  equation  (2-59).  Solving  for  xQ  n+^  imaginary  solutions 
appear  in  the  range  mentioned  above.  The  Carleman  approximation  method 
fails  to  predict  the  region  of  t where  random  oscillations  occur. 
Figure  2-2  compares  the  exact  solution  to  the  solutions  obtained  by  the 
proposed  new  method  over  one  period  after  the  transients  have  died 
out.  For  small  values  of  the  period  t near-singularity  of  the  matrix  to 
be  inverted  occurred,  so  numerical  solution  with  maximum  pivot-strategy 
was  employed  (see  Appendix  C).  Figure  2-2  clearly  shows  that  the 
approximation  to  the  exact  solution  by  the  new  method  is  excellent. 
Using  a third  order  Carleman  approximation  is  superior  to  second  order 
as  expected.  Table  2-1  shows  the  results  for  the  first  three  orders  of 
approximation.  Figure  2-3  shows  how  the  performance  measure  J varies 
with  the  period  t.  The  agreement  between  the  exact  solution  (equation 
(2-59) > and  the  approximations  is  good  for  small  values  of  the 
period  t and  becomes  worse  when  t approaches  the  region  of  chaotic 
oscillations.  This  is  to  be  expected  because  the  approximations  fail  to 
predict  the  existence  of  these.  Agreement  is  achieved  again  for  values 
of  x beyond  the  region  of  chaotic  oscillations. 

The  other  important  result  that  can  be  seen  in  Figure  2-3  is  that 
the  performance  measure  declines  with  increasing  period  x.  This 
predicts  that  steady-state  operation  (continuous  harvesting)  is  superior 
to  periodic  harvesting.  If  continuous  harvesting  is  not  possible,  the 
time  intervals  between  harvests  should  be  kept  as  short  as  possible. 


26 


Figure  (2-1).  Poincare  Mapping  Procedure  for  the  Verhulst  model 
(a=l)  a)  X=2,  b)  x=5,  c)  x=8.  As  x increases,  the 
y=x-lime  shifts  relatively  to  the  left.  For  x=2 
and  x=S  the  absolute  value  of  the  slope  at  the 
point  of  intersection  is  less  than  1,  thus  guaran- 
teeing. a stable  x-periodic  solution.  For  x=5  it  is 
greater  than  1,  so  no  convergence  can  occur  [36]. 


x (time) 
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Figure  (2—2) . Effect  of  the  order  of  the  Carleman  approximation. 

After  all  the  transients  have  died  out,  the  population 
x is  plotted  versus  time  for  one  period.  The  exact 

solution  ( ) is  compared  to  the  second  ( ) and 

third  order  approximation  (x=l). 
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Table  2-1.  Results  for  the  State  Vector  and  the  Performance  Measured 
for  the  First  Three  Orders  of  Approximation 


1st  order 


2nd  order 


3rd  order 


x , = x = m - m(e  T-1 ) 1 

v p I O 


'i  - 7 In  - e'T)l  x0>1 


X „ = X 

0,1  0 


X = 

0,2  o 


-2x  2 -t 

-me  + m m e 

B^e"1  - 1)  B1 


2m 


m 


B1 (e~T  - 1)  B1 


■ T )('-e'T)X0,1  * (e 


■t  _ e 2t  l 

2 “ 2 Xo,2 


X = x = 

0,1  o 


-mB^e  2t  + mB^  + 3m^B^(e  2t  + 1 ) 
(e“T  - 1) 


2 *”  T Q __  q 

- B m e + 3nTB.(me  T + 1)  B.nr 

+ i— + _I 


B1B3 


B_ 


[2]  2m2B1B3  + 3m3B2(e'2T  + 1) 

0,2  ° B2B3(e'T  - 1) 


B^m2  + 3m3B2(me  T + 1)  B2m3 

BfB3  B1B3 
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Table  2-1.  continued 


where 

B1  = e”T(e~T  - 2m)  - 1 

B2  = [2e"T(e"T  - 1)](e“T  -m) 

B.  = e"3T  + 3m2e“T(e'T  - 1)  - 1 - ^e'T  + 1)(e"T  ~ m)(e"T  ~ 1 >e 

e'T(e"T  - 2m)  - 1 


1) 
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Effect  of  the  period  x on  the  performance  measure  J. 

The  exact  solution  ( ) is  compared  to  the  first 

(•••)»  second  ( ) and  third  order  approxima- 

tion. 


Figure  (2-3). 
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Conclusions 

The  new  method  proposed  allows  to  solve  for  optimal  impulsed  forc- 
ing analytically  by  employing  the  Carleman  linearization  procedure.  The 
quality  of  the  approximation  is  excellent,  depending  on  the  order  of  the 
Carleman  linearization  procedure  used.  This  method  is  limited  to  sys- 
tems or  regions  of  systems,  where  inflicting  disturbances  with  a period 
t lead  to  x-periodic  output. 


CHAPTER  III 

APPLICATIONS  IN  BIOCHEMICAL  ENGINEERING 

III.  1.  Periodic  Optimization  of  Continuous 
Microbial  Growth  Processes 


Introduction 

Very  little  work  in  biochemical  engineering  has  been  done  with 
continuous  reaction  systems.  And  not  surprisingly  even  less  work  has 
been  done  in  control  and  optimization  of  biochemical  reactors.  However, 
over  the  past  few  years  continuous  fermentations  have  dramatically 
increased  in  significance  and  hence  their  optimization  becomes  very 
important . 

As  mentioned  in  Chapter  I,  periodic  reactor  operation  is  often 
superior  to  steady-state  operation,  which  gives  great  significance  to 
the  problem  of  periodic  optimization  of  continuous  microbial  growth 
processes . 

The  objective  of  this  work  is  to  apply  the  ir-criterion  to  micro- 
bial growth  systems  in  a continuous  chemostat  culture.  In  particular 
two  different  problems  are  addressed.  The  first  consists  of  optimiza- 
tion of  the  total  biomass  productivity,  based  on  dynamic  models  of  the 
chemostat.  The  second  problem  to  be  addressed  is  the  maximization  of 
the  productivity  of  protein  based  on  a structured  model  due  to  Williams. 
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Periodic  Optimization  of  Biomass  Production 

The  proper  choice  of  the  performance  measure  is  perhaps  the  most 
important  decision  in  any  optimization  process.  For  processes  that  have 
as  an  objective  the  production  of  biomass,  the  total  biomass  producti- 
vity (in  grams  per  unit  time  per  unit  volume)  as  represented  by  the 
product  Dx,  where  D is  the  dilution  rate  and  x the  biomass  concentra- 
tion, should  be  included  in  the  performance  measure.  Since  the  limiting 
substrate  constitutes  the  major  running  cost,  a reasonable  performance 
measure  is 


J = Dx  - qD3°  (3-1) 

where  3°  is  the  feed  substrate  concentration,  and  q is  the  relative  cost 
of  the  substrate.  All  mathematical  models  of  the  chemostat  that  have 
appeared  in  the  literature,  whether  unstructured  [38-40]  or  structured 
[41-48],  predict  that  J is  maximum  at  a unique  value  of  the  dilution 
rate  D.  This  optimum  dilution  rate  lies  between  0 and  the  maximum 

specific  growth  rate  at  the  prevailing  culture  conditions.  The 
existence  of  an  optimal  dilution  rate  has  also  been  verified 
experimentally  [49-52]. 

The  question  arises  whether  periodic  variation  of  the  dilution  rate 
can  lead  to  enhanced  performance  when  compared  with  the  optimum  steady- 
state  operation.  This  possibility  is  investigated  in  the  current  work 
employing  both  structured  and  unstructured  models. 

The  simplest  possible  model  type,  an  unstructured  model,  has  inher- 
ent the  assumption  of  "balanced  growth"  [53].  Following  this  assump- 
tion, the  dynamics  of  the  chemostat  can  be  described  by  the  following 
general  two-dimensional  system  of  equations: 
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— = p(s ,x)x  - Dx 


(3-2) 


ds 

dt 


(3-3) 


where  x and  s are  the  concentrations  of  biomass  and  substrate  respec- 
tively, p is  the  specific  growth  rate  and  Y is  a yield  factor.  Most  of 
such  unstructured  models  (with  the  exception  of  a model  due  to  Contois 
[54])  assume  that  n is  a function  of  s only.  By  far  the  most  commonly 
used  model  is  that  due  to  Monod  [38].  It  assumes  the  following  func- 
tionality for  the  specific  growth  rate: 


where  ym  is  the  maximum  growth  rate  at  the  culture  temperature  and  pH, 
and  Ks  is  a constant . For  the  Monod  model  (equations  (3-2)  to  (3-4)) 
the  optimal  dilution  rate  for  steady-state  operation  is  given  by 


(3-4) 


(3-5) 
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Consider  now  the  periodic  optimization  problem.  It  is  assumed  that  D is 
varied  periodically  with  an  average  equal  to  DQ  The  ir-criterion  can 
be  used  to  test  the  desirability  of  such  operation.  In  this  case  equa- 
tion (2-1)  takes  the  form  of  equations  (3-2)  and  (3-3)  with  p(s)  given 
by  (3-10*  u=D  (the  control  variable),  x^=x,  X2=s  (the  state  variables) 
and  y=x  (i.e.  it  is  assumed  that  the  biomass  concentration  is  the  mea- 
sured variable).  The  performance  measure  is  given  by  equation  (3-1), 
and  an  equality  constraint  of  the  form  (2-21)  is  relevant  with 

v = DS°  - c (3-6) 


where  c is  a constant.  This  requires  that  the  total  average  feed  is 
fixed.  The  Hamiltonian  function  in  this  case  is 

H - DX  - ,D3°*  - DX)  . »2<-  . D(3°  - 3>) 

s s 

+ y(DS°  - c) 

(3-7) 

Using  equations  (3-2)  - (3*7)  with  q=0,  tt(u)  is  found  to  be 
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where  all  variables  are  evaluated  at  their  steady-state  values  for 
D=Dopt*  As  can  be  easily  seen  from  equation  (3-8),  ir(o>)  is  negative  for 
all  values  of  the  cycling  frequency  w.  This  means  that  for  a Monod-type 
model  steady-state  operation  is  always  superior  to  cycling,  regardless 
of  the  cycling  frequency. 

All  unstructured  models  predict  a response  to  step-changes  in 
operating  variables,  such  as  3°  and  D,  which  is  faster  than  experimen- 
tally observed  [55].  This  is  a result  of  the  inherent  assumption  of 
those  models  that  there  is  no  time-lag  between  changes  in  the  substrate 
level  and  adjustment  of  the  growth  rate  at  the  appropriate  level.  To 
relax  this  assumption  one  may  assume  that  the  specific  growth  rate  is  a 
function  not  only  of  the  present  substrate  level  but  also  of  previous 
levels  in  a weighted  manner.  In  this  case  the  growth  rate  is  taken  to 
be  a function  of  z instead  of  s,  where 

z = j_m  s (t )F(t-i )di  (3-9) 

F(t-x)  is  a memory  function  that  weighs  appropriately  previous  sub- 
strate levels.  If  one  chooses 

F(Y)  = ae  (3-10) 

then  it  can  be  shown  [ 3^ U that  the  dynamic  equations  for  the  chemostat 
take  the  form 


= p(z)x  - Dx 


(3-11) 
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ds 

dt 


u(s)x  . , o 

— Y — + D (S  - s) 


(3-12) 


(3-13) 


A similar  delay-model  was  introduced  by  Wang  and  Stephanopoulos  [56]. 
The  smaller  the  magnitude  of  the  parameter  a,  the  slower  is  the  response 
of  the  growth  rate  to  a change  in  the  level  of  the  substrate.  Defining 
dimensionless  variables 


the  dynamic  equations  take  the  form 


dx 


zx 


(3-14) 


Dx 


dt  K + z 


s 


ds 


sx 


♦ D ( 1 - s) 


(3-15) 


dt 


dz  V % 
— = a(s  - z) 

dt 


(3-16) 


The  performance  measure  takes  the  form 


J = Dx  - D(1  - S) 


(3-17) 


where  g = 1 
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The  optimal  dimensionless  dilution  rate  is 


K 

s 

A 

K +8 
s 


The  Hamiltonian  in  this  case  becomes 


(3-18) 


H = Dx  - D(1  - s)  + A ( - Dx)  + A (-  ^ + D(1  - a)) 

K + z K + s 

3 s 

A A A A A 

+ A^a(s  - z)  + n(D  - c)  (3-19) 

The  expression  for  v(w),  which  is  again  a scalar  function,  is  given  in 
Appendix  D.  It  is  a function  of  only  three  variables  K , a and  g. 

To  test  the  effect  of  the  delay  a and  the  cost  factor  6,  it  was 
evaluated  numerically  and  plotted  as  a function  of  w.  When  q=0  (imply- 
ing B=1 ) , the  effect  of  delay  alone  may  be  seen  in  Figure  3-1.  It  is 
seen  that  nr(w)  becomes  positive  for  a range  of  cycling  frequencies 
u),  resulting  in  enhanced  performance.  It  was  also  found  that  the  great- 
er the  time  delay  (i.e.  the  smaller  a),  the  higher  is  the  maximum 
of  ir(ai)  and  the  lower  the  value  of  ui  at  which  it  occurs.  Large  values 
of  a result  In  an  all-negative  ir(w)  in  agreement  with  the  Monod  model 
without  delay  (equation  3-8).  Figure  3-2  shows  the  effect  of  the  dimen- 
sionless cost  factor  &.  As  8 increases,  the  maximum  of  ir( to)  decreases 
and  occurs  at  higher  cycling  frequencies.  Therefore,  higher  nutrient 
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tt(w)  28.0 


Figure  (3-1).  Tttto)  versus  to  for  Monod— model  with  delay  as  given 
by  equation  (3-10)  and  cost  factor  (dimensionless) 
Model  parameters : Kg  =0.9,  3=  1.,  a=  0.1,0. 5 
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Figure  (3—2).  tt(w)  versus  U)  for  Monod-raodel  with  delay  as  given 
by  equation  (3-10)  and  cost  factor  (dimensionless). 
Model  parameters:  £,=0.9,  a=l. , 6=0.1,  0.9. 


41 


costs  improve  the  effectiveness  of  cyclic  operation.  Notice  also  from 
equation  (3-18)  that  whenever  q is  greater  than  Y,  the  optimum  dilution 
rate  is  zero,  implying  that  in  this  case  a batch  operation  is  of 
choice.  Finally,  it  should  be  noted  that  for  very  large  cycling  fre- 
quencies, tt ( u> ) goes  to  zero,  implying  that  the  cyclic  operation  results 
in  an  average  performance  measure  equal  to  that  of  the  optimum  steady- 
state  operation.  This  is  expected  because  for  very  fast  cycling,  the 
culture  does  not  have  the  time  to  respond  to  periodic  variations,  so  it 
reacts  to  the  average  D,  which  is  D0p^ , and  the  performance  measure  J 
approaches  the  optimal  steady-state  value  J°  as  go  ■* 

If  in  place  of  equation  (3-10)  the  memory  function 

F(Y)  = S(y  - t)  (3-20) 

is  chosen  (6  being  the  Dirac  delta  function),  a somewhat  different 
behavior  is  observed.  In  this  case  growth  depends  on  the  substrate 
level  at  a particular  instant  in  the  past.  The  state  equations  are 
equation  (3-11)  and  (3-12)  with 

z = 3(t  - t)  (3-21) 

This  model  is  a two-dimensional  system  of  delay-differential  equa- 
tions. For  such  systems,  a ir-criterion  was  derived  by  Sincik  and  Bailey 
[23].  Their  results  were  used  to  find  an  expression  for  irU)  which  is 
given  in  Appendix  E.  For  this  delay  model,  more  than  one  frequency 

range  is  found  for  which  ir(o>)  > 0.  The  situation  is  exemplified  in 
Figure  3~3.  This  behavior  is  in  agreement  with  some  experimental 


findings  [15]. 
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Figure  (3-3).  tt(oj)  versus  0)  for  Monod  model  with  delay  as  given 

by  equation  (3-20).  Model  parameters:  y =0.7,  Y=0.5, 
Ks=22.,  S =30.,  t=  0.3  m 
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To  verify  the  predictions  of  the  ir- criterion , numerical  simulations 
were  undertaken.  Square-wave  variations  (bang-bang  switches)  around  the 
optimal  dilution  rate  were  employed,  rather  than  sinusoidal  variations, 
since  such  type  of  cycling  has  previously  proven  superior  and  easier  to 
implement  experimentally.  Figure  3~4a  is  a plot  of  the  ratio  of  perfor- 
mance measures  ^cyciing/Jsteady-state ^ versus  cycling  frequency  for  the 
system  described  by  equations  (3-1 4)- (3-1 6)  (with  a - 0.1,  B - 0.65 
and  Kg  = 0.5).  The  amplitude  of  the  square-wave  was  chosen  to  be  0.5 
Dopt  ’ For  comparison  purposes  a plot  of  ir(u>)  is  also  shown  (Figure  3- 
4b).  The  simulations  predict  remarkable  improvement  in  process  perfor- 
mance. Moreover,  the  optimal  cycling  frequency  as  obtained  by  intergra- 
tion  is  seen  to  be  close  to  that  predicted  by  the  ir-criterion. 

Periodic  Optimization  of  Protein  Production 

In  the  previous  section  time-delay  models  were  used  to  describe  the 
dynamics  of  the  ehemostat.  An  alternative  type  of  structured  modeling 
which  realizes  that  under  dynamic  conditions  balanced  growth  is  not 
valid  could  be  used  instead.  Such  structured  modeling  has  been  found, 
in  general,  superior  in  describing  the  dynamics  of  the  ehemostat.  The 
simplest  structured  model  is  that  due  to  Williams  [41]. 

The  fundamental  assumption  of  the  Williams  model  is  that  the  cell 
consists  of  two  basic  components,  one  synthetic  and  one  structural/gene- 
tic. The  synthetic  portion  includes  the  sum  total  of  small  metabolites 
in  the  cell  and  the  ribosomes  (RNA),  whereas  the  greater  bulk  of  the 
structural/genetic  portion  is  protein  and  DNA.  Although  the  exact 
physical  significance  of  the  variables  in  the  Williams-model  has  not 
been  established,  the  predictions  of  this  model  are  still  worth 
considering.  Williams's  model  in  dimensionless  form  may  be  written  as 
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Figure  (3-4).  Predictions  for  the  optimal  cycling  period 

a) J(to)  versus  to  for  Monod— model  with  delay  as  given 

by  equation  (3-10)  and  cost  factor  (dimensionless). 
Model  parameters:  K =0.5,  =0.1,  =.65,  D =o.34. 

b) TT(a))  versus  uj.  Model  parameters  as  in  a).  °^t 
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dx 


(3-22) 


df 


x^2  - Dx 


dx 


2 


(3-23) 


dt 


x^2  + (1  - x2)D 


Dx 


(3-24) 


dt 


a 


3 


where  x^  is  the  dimensionless  biomass  concentration,  x2  is  dimensionless 
substrate  concentration  and  x^  is  the  dimensionless  concentration  of  the 
structural/genetic  portion.  In  particular 


where  K1  and  K2  are  kinetic  constants  of  the  original  Williams  model. 

The  structural/genetic  portion  consists  mainly  of  protein,  the  DM 
content  being  relatively  small.  Thus  we  can  loosely  assume  that  x^  is 
the  dimensionless  concentration  of  protein  in  the  culture. 

We  are  interested  in  operating  the  chemostat  in  a manner  which  is 
optimum  for  the  production  of  protein.  That  is  in  this  case  the  perfor- 


mance measure  is 
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J = 


(3-25) 


The  optimum  dilution  rate  for  steady-state  operation  in  this  case  is 
calculated  to  be 


1 

2(1  + a ) 


(3-26) 


The  optimization  is  again  subject  to  the  equality  constraint  (3-6). 

Taking  u=D,  >c  = (x^  , , x^)^"  y*x^,  the  Hamiltonian  in  this  case 

becomes 


H = Dx3  + A1(x1x2  - Dx^  + ^2 (_Xl X2  + 0"X2)D)  + 
x (x  - x ) 

*3(  = Dxo)  + m(D  - c)  (3-27) 

a 3 

tt(o))  is  again  a scalar  and  the  appropriate  expression  is  given  in  Appen- 
dix F.  A graphic  illustration  is  given  in  Figure  3-5.  From  this  figure 
it  can  be  seen  that  cycling  results  in  improved  performance.  It  should 
be  noted  that  as  a increases,  the  optimum  forcing  frequency  decreases. 
Numerical  simulations  have  verified  the  global  validity  of  these  re- 
sults . 


Concl us ions 

The  n-criterion  is  a very  useful  method  for  determining  whether 
periodic  operation  can  lead  to  superior  bioreactor  performance.  En- 
hanced productivity  of  biomass  results  when  a bioreactor  is  operated 
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Figure  (3—5).  'n'(to)  versus  co  for  Williams'  model. 

5=0.6,  1.1 


Model  parameters 
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with  a periodically  varying  dilution  rate.  This  is  a consequence  of  the 
fact  that  there  is  a time  lag  in  the  culture  adaptation  to  a new  growth 
level  whenever  a disturbance  is  introduced.  Finally  structured  models 
such  as  the  Williams'  model  that  relax  the  assumption  of  "balanced 
growth,"  predict  an  increase  in  the  average  protein  productivity  upon 
periodic  variation  of  the  dilution  rate. 

Chapter  III.  2.  The  Effect  of  Periodic  Forcing  on 
Saccharomyces  cerevisiae  in  Continuous  Culture 


Introduction 

As  mentioned  in  Chapter  1,  continuous  operation  of  bioreactors  has 
only  become  important  over  the  past  few  years.  Thus,  although  a sub- 
stantial amount  of  work  has  been  undertaken  for  the  optimization  of 
batch  processes  [1-6]  only  very  little  has  been  done  for  the  optimiza- 
tion of  continuous  bioreactors  [7],  even  less  so  in  the  area  of  periodic 
reactor  operation. 

A few  workers  have  conducted  experiments  to  observe  the  behavior  of 
chemostat  cultures  when  subjected  to  cyclic  operation.  The  major  refer- 
ences are  summarized  in  Table  3~1  [15-20,  57,  58].  It  is  hard  to  draw 
general  conclusions  from  those  first  attempts  to  assess  the  effect  of 
cycling  experimentally.  In  many  cases,  complex  (e.g.  undefined)  growth 
media  like  molasses  were  used  or  some  quantities,  such  as  yield  factors 
have  either  been  not  well  defined  or  defined  differently.  In  these 
experiments  the  limiting  substrate  concentration  or  the  dilution  rate 
served  as  the  manipulated  process  variable. 

Picket,  Bazin  and  Topiwala  [15]  (using  S°,  the  feed  limiting  sub- 
strate concentration  as  control  variable)  found  that  at  some  frequencies 
periodic  reactor  operation  was  superior  to  steady-state  operation. 


Table  3-1.  Summary  of  Major  References  on  Experiments  Involving 

Periodic  Reactor  Operation 
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Moreover,  the  observed  improvement  in  these  frequency-ranges  was  found 
to  increase  with  increasing  cycling  amplitudes. 

It  has  been  known  for  some  time  that  the  concentration  of  some 
important  metabolites  inside  the  cell  is  not  constant  during  a dynamic 
process.  It  has  also  been  observed  that  the  average  macromolecular 
composition  of  the  cell  changes  under  cyclic  conditions  [17],  So  it  was 
found  that  the  percent  composition  of  protein  and  RNA  (essential  for 
growth)  is  higher  when  periodic  operation  is  employed.  Especially  in 
the  case  where  protein  is  an  economically  important  factor  this  result 
is  of  considerable  significance. 

The  purpose  of  this  study  is  to  observe  the  behavior  of  S.  cerevi- 
siae  in  continuous  culture  undergoing  square  wave  forcing  in  the  dilu- 
tion rate.  After  a steady-state  model  has  been  obtained  from  steady- 
state  data,  a dynamic  model  is  found  which  succeeds  in  describing  the 
response  of  the  system  to  step  changes  in  the  dilution  rate.  Subse- 
quently experiments  are  performed  to  observe  the  behavior  of  the  culture 
under  periodic  conditions. 

In  the  next  section  the  experimental  apparatus,  materials  and 
methods  will  be  described.  Thereafter  the  steady-state  as  well  as  the 
dynamic  response  experiments  are  described  and  it  is  shown  how  a mathe- 
matical model  was  derived  from  the  obtained  data.  The  last  part  of 
Chapter  3 is  concerned  with  the  periodic  optimization  of  the  bioreactor, 
in  theory  and  in  experiment. 

Equipment 

A schematic  representation  of  the  experimental  system  is  given  in 
Figure  3~6. 


Gas 
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Figure  (3-6).  A schematic  representation  of  the  experimental  system 
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The  Bioengineering  KLF  2000  fermentor  with  3 liter  glass  cylinder 
serves  as  chemostat.  The  average  reactor  volume  is  2.5  It.  A DC  motor 
driven  from  below  rotates  axially  two  flat  blade  agitators.  Mixing  is 
enhanced  by  steam  breakers  inside  the  vessel.  The  agitation  system  is 
connected  with  the  control  panel  for  speed  measurement  and  setting.  The 
rotation  speed  was  kept  at  500  rpm.  The  top  of  the  vessel  has  access 
ports  for  pH  measurement  and  control,  DO-electrode , foam  sensor  and 
control „ aeration,  addition  of  the  growth  medium  and  an  additional  port 
for  re-entry  of  a recycle  loop  for  the  turbidity  measurement . The 
bottom  plate  of  the  fermentor  has  two  ports,  one  for  harvesting  and  one 
for  the  recycle  loop. 

The'  temperature  is  controlled  via  a temperature  controller  located 
in  the  control  panel.  A pt  100  temperature  sensor  monitors  the  fermen- 
tation temperature  continuously.  An  800  W heating  finger  keeps  the 
temperature  at  30°C,  the  optimal  growth  temperature  for  yeast. 

The  pH  controller  also  located  at  the  control  panel  is  quasi 
steady-state  three-point  eliminator.  A silver/silver  chloride  electrode 
i3  inserted  through  the  top  of  the  vessel  and  monitors  the  pH  contin- 
uously. The  controller  keeps  the  pH  inside  the  fermentor  at  4.0  through 
actuating  two  peristaltic  pumps  that  deliver  either  NaOH  or  HC1  as 
needed . 

The  available  pumps  are  constant  flow  rate  pumps.  In  continuous 
steady-state  operation  of  the  harvesting  pump  is  set  for  a certain 
flowrate  and  not  controlled,  the  feed  pump  is  controlled  by  the  load- 
cell control  system  (on/ off). 

A turbidity  measurement  device  is  employed  for  the  continuous 
measurement  of  the  optical  density,  from  which  the  concentration  of 
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biomass  can  be  obtained.  The  TMK  turbidity  monitor  by  bioengineering  is 
equipped  with  a flow-through  cuvette  through  which  the  recycle  stream  is 
guided. 

A continuous  glucose  analyzer  from  Analytical  Research  was  used  for 
continuous  glucose  measurement.  Periodically  a small  amount  of  the 
medium  is  sampled  and  led  to  the  analyzer.  At  the  end  of  an  electro- 
chemical sensor  rests  an  immobilized  enzyme  membrane.  The  enzyme  cata- 
lyzes the  reaction  between  the  sugar  and  oxygen  to  form  H202  which  is 
detected  by  the  sensor. 

A PDP  11/23  microcomputer  equipped  with  a 10  M hard  disk  drive,  a 
graphics  terminal  and  an  LA  100  printer  together  with  Data  Translation 
AD/DA  converters  comprises  a complete  data  acquisition,  manipulation  and 
agitation  system,  to  be  used  on-line  for  direct  digital  control.  For 
the  case  of  periodic  squarewave  forcing  the  computer  has  been  programmed 
to  control  the  harvest  pump  in  an  on/ off  manner.  The  control  panel 
adjusts  the  feed  pump  to  keep  the  reactor  volume  constant. 

Materials  and  Methods 

S.  cerevisiae  was  obtained  from  the  Department  of  Microbiology  at 
the  University  of  Florida.  Stock  cultures  were  maintained  on  nutrient 
broth  agar  at  8°C.  The  composition  of  the  feed  medium  for  continuous 
operation  is  listed  in  Table  3-2.  To  avoid  the  possibility  of  a change 
in  the  feed  composition  through  a chemical  reaction,  batches  I and  II 
were  autoclaved  separately  and  mixed  under  UV-light.  For  the  start-up 
proceedings  batch  II  was  sterilized  in  the  fermenter  at  1 21 °C  for  30 
minutes,  batch  I was  autoclaved  and  used  as  a seed  flask  for  the  reac- 


tor . 
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Table  3.2.  Composition  of  the  Growth  Medium 


Substance 

Concentration 

Glucose 

2 g/l 

I 

Yeast  Extract 

.2  g/l 

Mg30i4*7H20 

.4  g/l 

II 

Potassium  Phosphate  monobasic 

5 g/l 

Ammonium  Sulfate 

2 g/l 

Dry  weights  of  cell  mass  were  obtained  by  heat-shocking  50  ml  of 
sample  for  6 minutes  to  inhibit  further  growth,  followed  by  centrifuga- 
tion (30  minutes  at  1200  rpm).  The  liquid  was  decanted  and  the  cells 
washed  with  distilled  water.  After  suspending  them  in  2 ml  of  distilled 
water  they  were  poured  in  open  petri-dishes  and  the  water  was  gently 
evaporated. 

Establishment  of  a Steady-State  Model 

The  material  balances  for  the  chemostat  can  be  described  by  the 
following  two-dimensional  system  of  equations  (unstructured  model) 

dx 

— = p(s)x  - Dx  (3-28) 

= - Y p(s)x  + D(S°-s) 


ds 

dt 


(3-29) 
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where  x and  s are  the  concentrations  of  biomass  and  substrate, 
respectively,  y is  the  specific  growth  rate  and  Y the  yield  factor 


grams  biomass  produced 
grams  glucose  used 


(3-30) 


The  unstructured  model  due  to  Monod  [38]  assumes  the  following  function- 
ality for  the  specific  growth  rate 


being  the  maximum  growth  rate  at  the  given  culture  conditions  and  K 

s 

being  a constant.  The  model  used  to  describe  the  steady-state  behavior 
of  the  system  therefore  is 


dx  Mm 


dt  K + s 
s 


x - Dx 


(3-3D 


ds  1 V 


dt  Y K + s 
s 


— X + D(3  -s) 


(3-32) 


At  steady  state  equation  (3-31)  can  be  rewritten 


as 


(3-33) 


This  functionality  can  be  used  to  determine  the  constants  y and  K by 

171  S 

plotting  1 /D  versus  1/s  (also  known  as  Lineweaver-Burk  plot).  Table  3-3 
summarizes  the  steady-state  measurements  of  the  system. 


56 


Table  3-3.  Summary  of  the  Steady-State  Behavior 


Dilution  Rate 
[hr-1 ] 

Biomass  Concentration 
[g/L] 

Substrate  Concentration 
[g/L] 

0.100 

.275 

.460 

-=J- 

O 

• 

o 

.290 

.403 

0.116 

.268 

.530 

0.154 

.173 

.802 

0.200 

— 

1 .196 

0.209 

.133 

1 .181 

Figure  3~7  shows  the  Lineweaver-Burk  plot  obtained  from  steady-state 
data.  Given  the  culture  conditions  described  in  the  previous  sections 
and  using  linear  regression  the  constants  were  found  to  be 


u = 0.5574 
m 

K = 2.0535 

3 

Y =0.175 

The  Response  of  the  System  to  a Step  Change  in  the  Dilution  Rate 

The  Monod-model , although  it  succeeds  in  describing  the  chemostat 
at  steady-state  in  most  cases,  is  known  to  fail  to  predict  dynamic 
responses  accurately.  Figure  3-8  proves  this  to  be  true  for  the  system 
at  hand.  The  response  of  the  biomass  in  graras/liter  is  followed  over 
time  after  a step  change  in  the  dilution  rate  from  0.2  to  0.1  has  been 
performed.  As  can  be  seen  clearly,  the  Monod-model  is  not  only  too 
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Figure  (3-7).  Lineweaver  - Burk  plot  for  the  determination  of  the  constants 
K and  p in  the  Monod-model 

g ' m 
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TIME  (min) 


•Kigure  (3-8).  The  response  of  the  biomass  to  a step-change  in 
the  dilution  rate  from  0.2  to  o.l;  the 
predictions  of  the  Monod-model 
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sluggish  to  portray  the  fast  initial  response  of  the  culture,  it  also 
fails  to  predict  the  overshoot  and  subsequent  return  to  the  steady-state 
value  corresponding  to  the  new  dilution  rate.  This  is  due  to  the  inher- 
ent assumption  of  this  model,  that  there  does  not  exist  a time-lag 
between  a change  in  the  substrate  concentration  and  adaptation  of  the 
culture  for  growth  at  the  new  substrate  level. 

To  relax  this  assumption,  a time-delay  was  introduced  in  the  model 
(see  Chapter  III.1  .).  Equation  (3-9)  was  used  and  the  best  fit  to  the 
experimental  data  was  found  for  the  case  of  a pure  time-delay,  i.e. 
choosing  equation  (3~20)  as  memory  function.  An  equivalent  way  to 
formulate  the  pure  time-delay  in  a form  which  is  convenient  for  numeri- 
cal simulations  is  (for  high  enough  k [59]) 


dx 

dt 


RZ  X 

m 1 

k + z„ 
s 1 


Dx 


ds 

dt 


1 vx  + 

y k3  + 3 + 


D(S°-s) 


) 


) 


(3-34) 


(3-35) 


(3-36) 


5zk 

dT"  = a(s  V 

a being  an  adaptability  parameter.  It  was  found  that  k = 10  is  suffi- 
ciently large.  Given  the  Monod-parameters  obtained  in  the  previous 
section,  it  was  attempted  to  find  the  value  of  the  adaptability 
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parameter  a which  would  give  the  best  fit  to  the  experimental  data,  as 
judged  by  the  "least-square-error If  one  tries  to  obtain  the  best  fit 
to  all  the  data  until  the  transients  have  virtually  subsided  (1,020 
minutes),  the  optimal  a is  found  to  be  3.  As  seen  in  Figure  3“9  the  fit 
is  not  in  very  good  quantitative  agreement  with  the  experimental  step 
response,  especially  in  describing  the  initial  steep  increase,  but  such 
a delay  model  is  able  to  predict  a faster  rise  in  biomass  than  the 
Monod-model.  This  is  a result  of  the  specific  growth  rate  continuing  to 
be  high  for  a certain  period  of  time,  following  a step  down  in  the 
dilution  rate.  On  the  other  hand,  if  the  objective  is  to  describe  the 
initial  steep  response  as  accurately  as  possible,  a much  lower  value 
for  a,  i.e.  1,  should  be  chosen.  But,  as  Figure  3“10  indicates,  using 
this  high  value  for  the  time-delay  predicts  an  overshoot  in  the  biomass 
concentration,  which  is  unacceptable. 

From  the  physiological  point  of  view,  it  is  well  established  that  a 
culture  enters  a so-called  "lag-phase"  when  it  is  subjected  to  different 
environmental  conditions.  It  utilizes  this  time  to  adapt  its  metabolism 
to  the  new  environment.  After  this  adaptation  has  taken  place  it 
responds  much  more  actively  to  the  environment.  Incorporating  this 
effect  in  the  model  allows  to  describe  the  step  response  successfully. 
Using  the  model  as  described  by  equations  (3-34)  to  (3-36)  and 

a = 1 . for  t < 3 hours 

a = 4.  for  t 2 3 hours  * * 

gives  the  result  as  seen  in  Figure  3— 11. 

Mow  consider  a step  change  in  the  dilution  rate  in  the  opposite 
direction,  from  0.1  to  0.2.  Figure  3-12  shows  once  more  that  the  Monod 
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Figure  (3  9).  The  response  of  the  biomass  to  a step— change  in 
the  dilution  rate  from  0.2  to  0.1;  the 
predictions  of  a pure  time-delay  model  ( ct=3.) 
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TIME  (min) 


Figure  (3-10) . The  response  of  the  biomass  to  a step-change 
in  the  dilution  rate  from  0.2  to  0.1;  the 
predictions  of  a pure  time-delay  model  ( a=l.) 


CELL  CONCENTRATION  (g tf) 
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figure  (3-11).  The  response  of  the  biomass  to  a step-change  in 
the  dilution  rate  from  0.2  to  0.1  ; the 
predictions  of  the  model  which  incorporates 
a change  in  the  physiological  state  of  the 
culture 
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TIME  (min) 


Figure  (3-12) . The  response  of  the  biomass  to  a step-change 
in  the  dilution  rate  from  0.1  to  0.2;  the 
predictions  of  the  Monod-model 
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model  describes  a response  which  is  much  too  slow.  But  in  this  case 
(Figure  3~13)  even  the  pure  time-delay  model,  which  incorporates  a shift 
in  the  physiological  state  of  the  culture,  although  predicting  a faster 
response,  does  not  succeed  in  accurately  portraying  this  step  re- 
sponse. After  analyzing  the  data  it  was  concluded,  that,  after  the 
dilution  rate  was  stepped  up  from  0.1  to  0.2,  all  growth  has  stalled. 
As  Figure  3-1^  proves,  the  drop  in  the  biomass-concentration  corresponds 
exactly  to  the  one  expected  for  zero  growth  rate,  leaving  only  the 
hydraulic  effect  of  the  culture  being  washed  out  of  the  fermenter. 

It  might  seem  surprising  at  first  that  growth  should  be  arrested 
following  a step  up  in  the  dilution  rate,  but  this  state  of  shock  can  be 
compared  to  the  one  where  an  inoculum  is  placed  into  new  medium.  Again, 
the  culture  reacts  by  entering  a lag-phase  to  adapt  to  the  new  environ- 
ment. As  the  previous  figures  indicate,  growth  starts  to  reoccur  after 
a certain  time.  And,  not  surprisingly,  this  time  is  about  three  hours, 
thus  agreeing  with  the  length  of  the  lag-phase  for  a step  change  down  in 
dilution  rate.  The  length  of  the  lag  phase  very  likely  is  a function  of 
the  magnitude  of  the  step  change,  being  virtually  zero  for  small  steps 
and  higher  for  larger  ones. 

The  Response  of  the  System  to  Periodic  Variation  of  the  Dilution  Rate 

Starting  from  a steady-state  corresponding  to  a dilution  rate  of 

0. 1,  the  system  was  subjected  to  square  wave  variations  in  the  dilution 
rate  of  different  periods.  The  amplitude  was  twice  this  dilution  rate; 

1. e.  variations  from  D = 0.2  to  D = 0.0  (batch  operation)  were  implemen- 
ted, so  that  the  total  amount  of  substrate  over  a period  equals  the  one 
at  steady-state  operation  with  D = 0.1. 


CELL  CONCENTRATION  (g/i) 


66 


Figure  (3-13).  The  response  of  the  biomass  to  a step-change 
in  the  dilutionrate  from  0,1  to  0.2;  the 
predictions  of  the  model,  which  incorporates 
a change  in  the  physiological  state  of  the 
culture 
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Figure  (3-14).  The  response  of  the  biomass  to  a step-change 
in  the  dilution  rate  from  0.1  to  0.2;  the 
response  as  it  would  be  observed  for  zero 
growth  rate  (washout) 


68 


From  the  step-response  experiments,  which  indicated  a time-lag  of 
about  three  hours,  it  was  concluded  that  a cycling  period  in  the 
vicinity  of  three  hours  might  influence  the  performance  of  the  syste 
positively.  Applying  the  model  established  in  the  previous  section, 
qualitative  agreement  is  found.  Cycling  with  a period  of  two  hours  is 
expected  to  yield  lower  productivities  than  steady-state  operation, 
cycling  with  a period  of  three  hours  higher  ones.  This  modeling  is  not 
applicable  to  the  four-hour  period  because  the  model  does  not  account 
for  a possible  entry  into  the  stationary  phase  for  high  periods. 

Figure  3~15  presents  the  results  (after  all  transients  have  sub- 
sided). While  cycling  with  a period  of  either  two  or  four  hours  does 
not  produce  as  much  biomass  as  steady-state  operation,  cycling  with  a 
period  of  three  hours  results  in  significantly  improved  reactor  perfor- 
mance. Table  lists  the  concentrations  during  the  cycle  and  the 

productivities . 


Table  3-4.  Biomass  Concentrations  and  Productivities 
During  Cyclic  Reactor  Operation 


T 

Concentration  [g/fc] 
low  high 

Productivity  [g /hr] 

2 

0.15 

0.22 

0.0518 

3 

0.29 

0.40 

0.0952 

4 

0.09 

0.13 

0.0305 

steady  state 
(D  = 0.1) 

0 

.275 

0.0756 

CELL  CONCENTRATION  (g H)  CELL  CONCENTRATION  (git)  CELL  CONCENTRATION  (git) 
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Figure  (3-15) . Biomass  concentration  versus  time  for 
different  cycling  periods 
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The  improvement  observed  at  a cycling  period  of  three  hours  is 
quite  remarkable.  During  the  half  of  the  cycle  when  the  dilution  rate 
is  equal  to  zero,  essentially  batch  growth,  the  culture  clearly  is  in  a 
state  of  exponential  growth  and  during  the  remainder,  when  the  dilution 
rate  is  high,  growth  has  stalled.  This  confirms  the  observations  made 
in  the  step  change  experiments.  The  magnitude  of  the  improvement  might 
be  explained  by  the  fact  that,  in  this  case,  a high  (exponential)  growth 
rate  coincides  with  a high  concentration  of  biomass  at  a time  where 
there  is  ample  substrate  present. 

At  a cycling  period  of  four  hours,  as  Figure  3-15  shows,  the  cul- 
ture is  not  in  a state  of  exponential  growth  during  the  half  of  the 
cycle  when  batch  operation  prevails.  The  leveling-off  of  the  slope 
indicates  entry  into  the  stationary  phase,  while  once  again,  no  growth 
is  found  for  the  rest  of  the  cycle.  The  culture  thus  switches  between 
nearly  stationary  phase  and  lag-phase,  not  optimal  conditions  to  operate 
under.  That  explains  the  very  low  productivity  at  this  cycling  fre- 
quency . 

The  reduced  productivity  at  a period  of  two  hours  might  be  caused 
by  the  fact  that  the  culture  is  subjected  to  a change  in  the  environmen- 
tal conditions  too  frequently  and  does  not  have  time  to  recover  during 
the  repeated  shocks. 

Conclusions 

Subjecting  the  culture  to  a sudden  step  up  in  dilution  rate  causes 
it  to  enter  a lag- phase  during  which  growth  is  temporarily  halted.  A 
step  down  in  dilution  rate  has  a similar  effect;  a lag-phase  of  about 
equal  length  is  found  during  which  the  growth  rate  remains  close  to  that 
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before  the  change.  Thereafter  the  culture  adapts  to  the  new  environ- 
ment . 

The  same  responses  are  observed  when  the  system  undergoes  square 
wave  variations  in  the  dilution  rate.  The  reactor  performance  was  found 
to  be  inferior  to  steady-state  operation  when  a cycling  period  of  two  or 
four  hours  was  used,  but  at  a cycling  period  of  three  hours  a remarkable 
increase  in  the  biomass  productivity  was  obtained. 


CHAPTER  IV 

APPLICATIONS  IN  CANCER  CHEMOTHERAPY 
AS  AN  EXAMPLE  FOR  DRUG  THERAPY 

Surgery  proves  ineffective  in  over  50?  of  all  cases  of  cancer  where 
the  neoplasms  are  nonlocal! zed  or  diffuse.  And  so  far  the  success  of 
chemotherapy  and  radiation  therapy  has  been  limited.  The  limiting 
factor  in  cancer  chemotherapy  (or  radiation  therapy  for  that  matter)  is 
the  toxicity  of  the  antineoplastic  agent  to  the  normal,  non-carcinogenic 
cell  population.  The  current  approach  to  determination  of  a good  treat- 
ment for  an  individual  in  clinical  practice  is  still  a trial  and  error 
process.  Different  treatment  centers  use  one  and  the  same  drug  in  as 
many  different  ways.  Berenbaum  [59]  proved  that  a drug  successful  in 
reducing  the  cancer  population  while  keeping  the  normal  cell  population 
above  a certain  level  can  easily  be  classified  as  ineffective  if  given 
the  wrong  way  (even  if  the  overall  amount  of  drug  stays  the  same). 

The  long  overdue  methodical  search  for  the  optimal  treatment  regi- 
men can  only  be  done  by  drawing  the  complete  picture  which  takes  into 
account  both  tumor  and  normal  cell  kinetics  and  drug-cell-interaction  as 
well  as  drug  resistance  and  pharmacokinetics.  Given  just  limited 
information  about  the  growth  history  of  the  tumor  and  given  the  pharma- 
cological properties  of  the  drug,  a mathematical  model  will  be  able  to 
predict  the  response  of  tumor  and  normal  cells  to  chemotherapy  and  thus 
allow  the  choice  of  the  proper  drug,  dosage  and  timing. 
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Mathematical  modeling  of  multicellular  organisms,  especially  as 
complex  ones  as  the  human  body,  is  extremely  difficult  because  many 
processes  at  the  cellular  level  like  drug-cell  interactions  and  mass 
transfer  limitations,  just  to  name  a few,  are  still  not  well  under- 
stood. For  the  lack  of  available  data  and  parameters,  the  model  has  to 
be  kept  as  simple  as  possible.  Emphasis  is  put  on  clinical  application 
rather  than  detail.  The  long  term  goal  is  to  make  use  of  the  kinetic 
differences  between  normal  and  malignant  cell  populations  in  order  to 
develop  chemotherapeutic  protocols  with  increased  selectivity  in  killing 
the  neoplastic  cells.  The  purpose  of  this  investigation  is  to  establish 
guidelines  for  optimal  treatment  using  the  limited  information  available 
and  to  show  how  optimal  control  theory  can  be  a powerful  tool,  leading 
to  new  approaches  in  chemotherapy. 

IV. 1.  A Novel  Approach  for  Determining  Optimal 
Treatment  Regimen  for  Cancer  Chemotherapy 


Introduction 

It  has  frequently  been  questioned  whether  continuous  or  periodic 
treatment  renders  the  highest  effectiveness.  The  problem  of  determining 
the  optimum  periodic  operation  for  a process  described  by  ordinary 
differential  equations  was  addressed  by  Bittanti  et  al . [22]  who  deve- 
loped the  so  called  ir-criterion  for  optimality.  As  explained  in  Chapter 
II. 1.,  this  criterion  states  a sufficient  condition  for  improvement  of 
the  performance  of  a dynamic  system  via  cycling  and  in  addition  provides 
an  estimate  of  the  optimum  cycling  frequencies  for  the  particular  system 


of  interest. 
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The  objective  now  is  to  investigate  the  consequences  of  the 
v-criterion  for  cancer  chemotherapy.  A procedure  will  be  given  for 
determining  for  a certain  patient  and  a given  drug  whether  continuous  or 
periodic  treatment  will  give  better  results  and  for  finding  the  optimal 
time  on/off  the  drug.  The  next  section  sets  the  mathematical  framework 
for  studying  the  posed  optimization  problem.  The  subsequent  sections 
present  the  results  of  applying  the  ir-criterion  to  this  physiological 
system. 

Establishment  of  a Suitable  Mathematical  Model  and  Performance  Measure 

As  mentioned  earlier , the  limiting  factor  in  cancer  chemotherapy  is 
the  toxicity  of  the  drug  to  the  "limiting  tissue,"  normal  homeostati- 
cally  controlled  body  cells.  Although  all  body  cells  are  affected  by 
the  drug,  there  is  one  "limiting  tissue"  which  in  most  cases  is  the  bone 
marrow  or  the  intestinal  cells.  Therefore  two  cell  populations  have  to 
be  modelled,  the  neoplastic  and  the  limiting  tissue  cells. 

The  following  cellular  material  balances  may  be  written  for  both 
cell  types: 

dx. 

— — = f . (x. ) - g. (c)x.  i = c or  n (4-1 ) 

dt  11 

A A 

where  xfi  is  the  normal  and  xq  the  cancerous  population,  c is  the  drug 
concentration  and  t the  time.  Several  models  have  been  proposed  to 
describe  the  cellular  growth  rate  functions  f^x  ) [59-65].  For  most 
human  tumors  the  Gorapertz  equation  provides  the  best  fit 
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dJlnx^  x1 

— x — = A - a Jin  — A , a ...  constants  (4-2) 

dt  x 0 

o 

It  predicts  fast  exponential  growth  for  a small  tumor  and  slowing  growth 
for  increasing  tumor  size.  Several  authors  [66,67]  successfully  mod- 
elled tumor  growth  using  this  equation. 

The  drug-cell-interaction  function  g^c)  has  been  modelled  either 
as  a linear  or  a Michaelis-Menten  type  expression.  It  is  a reasonable 
assumption  that  a cell  population  will  exhibit  a certain  saturation 
effect  towards  a drug,  so 


~ k.c 

g.  (c)  = k.  , K ...  constants  (4-3) 

K + c 1 

Swan  and  Vincent  [68]  used  an  equation  of  this  form  to  model  the  data  on 
ten  patients  who  had  been  treated  according  to  the  MCP-protocol  (Melpha- 
lan,  Cyclophosphamide,  Prednisone).  They  assumed  a single  fictitious 
drug  to  represent  the  effect  of  MCP.  The  parameters  obtained  from  their 
data  will  be  employed  in  this  work. 

Several  authors  [61,68,69]  solved  optimal  control  problems  for 
cancer  chemotherapy.  Their  result  was  a desired  plasma  concentration  of 
the  drug  versus  time  at  the  action  site.  But  these  results  are  of 
rather  limited  practical  use  because  they  cannot  tell  a clinician  how 
much  medication  to  administer  to  a patient  or  at  what  time  intervals. 

An  applicable  result  requires  the  regimen  to  be  not  in  terms  of 
plasma  drug  concentration  but  in  terms  of  actual  dosage.  This  in  turn 
calls  for  a relationship  which  describes  the  fate  of  the  drug  once  it 
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enters  the  body.  The  proper  pharmacokinetic  equation  for  the  case  of 
constant  infusion,  treating  the  body  as  one  compartment  is  [70] 

a 

do  **  a 

— = u - 6c  6 ...  constant  (4-4) 

dt 

A A 

where  c represents  the  plasma-concentration  and  u the  dosage  regimen, 
which  is  the  manipulated  variable  for  the  model  adopted  in  this  work. 
Making  the  following  assumptions: 


• Gompertz  growth  of  untreated  tumor 

• Saturation  effect  of  drug-cell  interaction 

• First  order  pharmacokinetics 

• One  compartment  model 

The  dimensionless  model  is  as  follows  (Appendix  G): 
dy  Y c 

(4-5) 
(4-6) 
(4-7) 


dt 

dy, 

dt 


= -<*y„  ~ 


n 1 + c 


= - y - 


Y C 
C 


c 1 + c 


normal  cell  population 


tumor  cell  population 


dc 

— = u - 6c 


drug  concentration 


The  steady  state  this  system  approaches  is  given  by 
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y, 


n 


o 


(4-8) 


o 


o 


(4-9) 


y, 


c 


0 


6 + u 


o 


o u 
c = — 


(4-10) 


As  can  be  seen  easily  from  equations  (4-8)  - (4-10)  in  connection  with 


postulated . 

"Steady  state"  in  the  strict  meaning  translates  into  continuous 
infusion  for  this  problem.  But  given  the  large  time  constants  of  this 
system  it  can  be  considered  to  be  at  steady  state  if  medication  is  given 
daily  as  opposed  to  treatment  cycles  of  e.g.  30  days  or  more.  This 
treatment  regimen  will  be  referred  to  as  "continuous  treatment." 

The  proper  choice  of  the  performance  measure  is  the  most  important 
decision  in  any  optimization  process.  For  this  physiological  system  the 
objective  is  twofold: 

1.  minimize  the  malignant  cell  population 

2.  maximize  the  benign  cell  population 


Appendix  C,  if  no  drug  is  given  (u°  = 0)  x 


In  the  case  of  multiple  objectives  they  can  be  linked  to  a single  per- 
formance measure  by  using  'weighing  factors.'  However,  the  choice  of 
these  factors  requires  a value  judgement  by  the  scientist,  which  is 
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often  difficult  to  make  ("How  important  is  it  to  save  normal  cells 
compared  to  killing  tumor  cells?").  Therefore  it  is  of  utmost 
importance  to  choose  a performance  measure  such  that  the  weighting 
parameter  has  a physical  meaning.  It  is  easy  to  see  that  exp  [y.  (t ) ] 
represents  the  ratio  of  the  number  of  cells  at  time  t of  treatment 
divided  by  the  number  of  cells  in  the  absence  of  treatment  (the  maximum 
attainable  number  of  tumor  cells).  A proper  performance  measure  is  then 


J = e n - qe  C (4-11) 

q being  a "cost  factor."  its  physical  meaning  is  the  percentage  of 
normal  cells  the  scientist  is  willing  to  sacrifice  in  order  to  kill  1$ 
of  the  cancer  population.  This  number  might  be  quite  high,  for  example 
in  cases  where  bone  marrow  transplants  are  possible  and  very  low  in 
other  cases. 

Optimization  of  Cancer  Chemotherapy 

Assume  that  the  behavior  of  benign  cells/malignant  cells  and  the 
drug  can  be  described  by  equations  (4-5)  - (4-7)  and  the  performance 
measure  is  given  by  equation  (4-8).  The  performance  measure  j is  a 
function  of  the  manipulated  variable,  the  dosage  rate  u.  Figure  4-1 
shows  the  dependence  of  J on  the  steady-state  value  of  u.  As  can  be 
seen  clearly  there  is  an  optimal  steady-state  dosage  uQ  For  the  case 
of  continuous  treatment  the  optimal  (dimensionless)  dosage  rate  is  given 


by 


PERFORMANCE  MEASURE 
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Figure  (4-1) . Performance  measure  J versus  rate  of  dosage  u.  Model 

parameters:  a=10,  y =401,  y =320.8,  6=66.7,  q=0.4 

c n n 
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u = on5i.nC 

opt  aY  - Y - alnE, 
c n 


(4-12) 


Now  consider  the  periodic  optimization  problem,  i.e.  whether  enhanced 
performance  is  possible  upon  periodic  variation  of  the  dosage.  It  is 
assumed  that  the  dosage  u is  varied  periodically  with  an  average  equal 
t0  uopt  as  shown  in  Figure  4-2.  The  methods  of  Chapter  2 can  be  used  to 
test  the  desirability  of  such  operations.  In  this  case,  equation  (2-1) 
takes  the  form  of  equations  (4-5)  - (4-7)  with  x.]  = xn  * x2  = yc  > x3  = 0 
(the  state  variables)  and  y_  = x_  (the  outputs  being  the  same  as  the  state 
variables).  The  state  variables  in  dimensionless  form,  for  steady  state 
operation  (continuous  treatment),  are  given  by  equations  (4-8)  to  (4-10) 
by  setting  u°  = uQpt.  The  performance  measure,  as  described  in  equation 
(2-19),  is  given  in  this  case  by  equation  (4-11).  Relevant  constraints 
are 

y(x)  = y(o)  (outputs  are  also  x-periodic)  (4-13) 

u(t)  = (there  is  a maximum  allowable  dosage)  (4-14) 

1 f T 

- JQ  u(t)dt  = u (total  amount  of  drug  over  a period  (4-15) 

is  fixed) 


This  gives  rise  to  the  Hamiltonian  function 
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Figure  (4-2).  Rate  of  dosage  u versus  period  T. 
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y, 


n 


y c 
n 


H = e 


n 1 + c^  + X2  ^ yo 


1 + e 


(4-16) 


+ A_  (u  - 6c)  + p (u  - u 
->  opt 


*i  and  p representing  Lagrange  multipliers  for  the  equality  and  inequa- 
lity constraints  respectively. 

Using  this  system  of  equations  tt(io)  is  calculated  to  be 


ir(w)  is  a function  of  all  five  dimensionless  model  parameters.  Depend- 
ing on  the  choice  of  those  parameters,  three  regions  can  be  seen 
(Figures  4-3,  4-4).  These  regions  can  be  characterized  as  follows: 

Region  1:  uQpt  negative 

Region  2:  m(u))  positive 

Region  3:  ir(ui)  negative 


tt(oj)  = 


(a:2  + 62)(1  + c°)3 


o 


(4-17) 


GAMMA  N/  GAMMA  C 
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ALPHA 


Figure  (4-3).  The  three  regions  for  q=0.5. 


GAMMA  N/  GAMMA  C 
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ALPHA 


Figure  (4-4).  The  three  regions  for  q=0.2 
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For  certain  values  of  the  five  parameters  uQpt , as  given  by 
equation  ( 4- 12),  is  negative.  Because  J is  a monotonic  function  of  u, 
the  optimal  dosage  for  this  case  is  uQpt  = 0.  This  simply  means  that 
the  drug  chosen  is  inappropriate  and  another  medication  has  to  be 
considered  whenever  the  parameters  give  rise  to  a point  in  Region  1. 

If  they  give  rise  to  a point  in  Region  2 (h(oj)  positive),  then  the 
model  predicts  that  some  form  of  periodic  treatment  will  be  superior  to 
continuous  treatment.  How  to  find  the  optimal  length  of  the  treatment 
period  x as  well  as  the  number  of  days  on  and  off  medication  will  be 
explained  later. 

For  Region  3»  it ( tu ) is  negative  which  means  that  continuous  treat- 
ment (drug  given  daily,  no  recovery  periods)  will  render  the  best  re- 
sults . 

To  test  the  effect  of  the  five  dimensionless  parameters,  tt ( to ) was 
evaluated  numerically  and  plotted  as  a function  of  <u.  These  results 
are  shown  in  Figures  4-5  to  4-9.  It  is  seen  that  tt(uO  becomes  positive 
for  a range  of  cycling  frequencies,  indicating  that  periodic  treatment 
using  these  treatment  frequencies  yields  better  results  than  continuous 
treatment . 

From  the  computation  of  xr ( to ) certain  important  qualitative  observa- 
tions may  be  made.  The  length  of  the  treatment  cycle  is  influenced  by 
the  growth  of  tumor  cells  relative  to  normal  cells  (the  parameter 
a),  and  by  the  tumor  cell  kill  that  the  medication  causes  relative  to 
the  kill  of  benign  cells  (Y  , Y ) . An  increasing  ratio  of  Y /Y  , which 
means  increasing  kill  of  normal  cells  compared  to  malignant  ones,  re- 
quires decreasing  cycling  frequencies  (w  .).  This  in  turn  means 

opt 

longer  treatment  cycles  are  necessary  to  give  benign  body  cells  time  to 


Pi  (omega) 
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Figure  (4-5).  The  effect  of  the  drug  instability  6 on  TT  ( to) . Model 

parameters:  a=10,  y =401,  y =320.8,  q=0.4 

c n 


Pi  (omega) 
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OMEGA 


Figure  (4-6).  The  effect  of  the  'cost  variable'  q on  tt  ( oo) . Model 

parameters:  a=10,  y =401,  y =320.8,  6=66.7 

c n 


Pi  (omega) 
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Figure  (4-7).  The  effect  of  the  growth  related  variable  a on  tt  ( to)  . 

Model  parameters:  Yc=401,  Yn=320.8,  6=66.7,  q=0.4 


Pi  (omega) 
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OMEGA 


Figure  (4-8).  The  effect  of  the  ratio  of  y /y  on  tt(oj).  Model 
parameters:  a=10,  yc=401,  <5=66?7 , q=0.4 


Pi  (omega) 
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Figure  (4-9).  The  effect  of  the  tumor  cell  kill  variable  y or  tt(w). 

Model  parameters:  a=10,  yn/yc=0.8,  6=66.7,  q=0.4 
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recover.  The  length  of  the  optimal  treatment  cycle  is  also  influenced 

by  the  parameters  a and  Y itself  but  not  in  a monotonic  fashion.  An 

c 

increasing  value  of  a,  which  means  increasing  growth  of  cancer  cells 
compared  to  normal  cells  first  causes  the  optimal  treatment  frequency  to 
decrease,  then  to  increase  again  (the  value  of  the  extremum  is  influ- 
enced by  the  other  parameters).  The  same  holds  true  for  the  parameter 
Y . The  parameters  which  do  not  influence  the  length  of  the  treatment 
cycle  are  the  drug  stability  (6)  and  the  "cost  parameter"  q.  It  might 
seem  surprising  at  first  that  the  drug  stability  influences  neither  the 
type  of  operation  (cyclic  or  continuous  as  can  be  seen  from  equation  (4- 
17))  nor  the  optimal  cycling  frequency,  but  it  strongly  influences  the 
average  amount  of  drug  optimally  administered  (uQpt).  This  indicates 
that  if  a drug  is  highly  unstable  it  should  be  given  at  a higher  dosage, 
if  permissible,  not  in  shorter  treatment  cycles  (for  this  set  of  growth 
and  kinetic  parameters).  It  is  fortunate  that  the  "cost  parameter"  q, 
which  requires  a value  judgement  by  the  physician  and  is  not  based  on 
physical  data,  does  not  influence  the  length  of  the  treatment  cycle. 

To  verify  the  predictions  of  the  ir-criterion  numerical  simulations 
(integrations  of  equations  (4-5)  to  (4-8)  and  (4-11)  were  undertaken. 
Square  wave  variations  (bang-bang  switches)  between  u = 2uQpt  and  u = 0 
(see  Figure  4-2)  were  employed,  rather  than  sinusoidal  variations,  since 
such  type  of  cycling  has  previously  proven  superior  and  is  easier  to 
implement.  Figure  4-10  shows  plots  of  the  performance  measure  and  the 
function-  ir(ui)  versus  the  cycling  frequency  u).  It  can  be  seen  that  the 
optimal  cycling  frequency  as  obtained  by  integration  (value  at  which  J 
is  maximum)  is  close  to  that  predicted  by  the  -ir-criterion  (value  at 


which  it  is  maximum). 


PERFORMANCE  MEASURE  pj  (omega) 
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Figure  (4-10).  Tr(u))  and  the  performance  measure  J(w)  versus  treat- 
ment frequency  U).  Model  parameters  as  in  Figure  4-1. 
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Although  the  ir-criterion  determines  the  optimal  cycling  frequency, 
it  does  not  give  any  predictions  about  the  optimal  wave  form,  i.e.  the 
optimal  fraction  of  the  treatment  cycle  drug  should  be  administered. 
Keeping  the  amount  of  drug  over  a period  constant,  the  variable  e is 
defined  as  that  fraction,  while  1-e  is  the  fraction  of  the  treatment 
cycle  for  which  the  patient  stays  off  medication  (Figure  4-11).  The 
numerical  simulations  show  that  there  exists  an  optimal  e,  as  shown  in 
Figure  4-12.  An  example  is  given  in  Appendix  H.  For  a tumor  of  the 
given  growth  characteristics  and  a drug  with  the  given  properties  it 
suggests  a treatment  cycle  of  35  days,  3 of  those  on  medication  and  32 
off,  as  well  as  the  optimal  rate  of  dosage.  In  some  cases  the  optimal 
fraction  e ^ cannot  be  implemented  because  the  dosage  rate  u,  which  is 
given  as  the  ratio  of  u /e  is  limited  by  umax,  the  maximum  allowable 
doses  according  to  toxicity  criteria.  So  there  is  a constraint  placed 
on  e given  by 


emin 


opt 


u 


max 


(4-18) 


If  the  calculations  yield  an  e . smaller  than  e . , then  the  fraction 

opt  min 

of  the  period  on  medication  has  to  be  set  at  e 

min 


Conclusions 

The  method  proposed  can  help  to  answer  the  fundamental  questions  of 
treatment  strategy  for  a certain  patient  and  the  drug  under  considera- 
tion: (1)  Is  the  drug  appropriate?  (2)  Will  continuous  or  cyclic 

administration  give  the  best  result?  (3)  How  much  shall  be  adminis- 
tered and  in  which  time  intervals?  The  quality  of  the  predictions 


DOSAGE  RATE 
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Figure  (4-11).  Different  wave  forms,  £=1,  1/2,  1/4. 
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Figure  (4-12).  The  performance  measure  J versus  e,  the  fraction  of 
the  treatment  cycle,  a patient  receives  medication. 
Model  parameters  as  in  Figure  4-4. 
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depends  solely  on  the  quality  of  the  available  data.  Although  the 
mathematical  model  used  in  this  paper  is  rather  rudimentary,  its  predic- 
tions are  based  on  the  actual  patient  data.  If  more  is  known  about  the 
growth  history  of  the  tumor  and  the  pharmacokinetics  of  the  drug  used,  a 
more  sophisticated  model  can  be  constructed.  The  more  complete  the 
available  data,  the  better  the  mathematical  model  and  the  better  will 
the  predictions  about  the  optimal  treatment  regimen. 

One  obvious  advantage  of  a "simplistic"  model  of  the  form  used  in 
this  work  is  the  ease  of  determination  of  the  model  parameters  from 
available  data.  In  addition,  the  mathematical  complexity  added  by  a 
more  sophisticated  model  can  make  the  determination  of  the  optimal  drug 
administration  significantly  less  tractable.  Calculating  optimal  treat- 
ment regimen  using  even  very  simplistic  models  represents  an  improvement 
compared  to  the  trial-and'-error  method  still  used  in  clinical  prac- 
tice. Summarizing  it  shall  be  stated  that  the  proposed  method  has 
proven  useful  in  determining  which  treatment  regimen  is  optimal  for  a 
certain  patient  and  a given  drug.  It  not  only  clearly  separates  the 
regions  in  which  continuous  or  periodic  treatment  yields  superior  re- 
sults, it  also  gives,  in  combination  with  numerical  simulations,  the 
optimal  length  of  the  treatment  cycle  as  well  as  the  optimal  time  on  or 
off  the  drug  and  the  required  dosage. 

IV. 2.  Periodic  Impulse  Forcing  in  Cancer  Chemotherapy- 
Determination  of  the  Optimal  Treatment  Regimen 


Introduction 

As  concluded  in  the  previous  Chapter  (IV. 1.),  optimal  control 
theory  in  general  and  the  ir-criterion  in  particular  can  make  a valuable 
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contribution  in  determining  how  promising  a certain  treatment  approach 
is.  The  appropriateness  of  a drug  (combination  of  drugs)  can  be 
established  as  well  as  the  form  of  treatment-continuous  or  periodic. 
The  ir-criterion,  as  explained  in  Section  II. 1,  was  applied  to  answer  the 
questions  stated  above. 

A point  of  particular  interest  is  how  the  body,  normal  and 
cancerous  cells,  reacts  to  a single  injection  of  drug  or  a series  of 
injections.  An  injection  of  a drug  into  the  body  represents  an 
instantaneous  addition  of  a substance  to  the  system  and  this  is  to  be 
described  mathematically  as  an  impulse.  Consequently  repeated 
injections  constitute  impulse  forcing.  Impulse  forcing  is  a very  common 
periodic  waveform  but  the  predictions  of  the  ir-criterion  are  less 
reliable  here,  mainly  for  two  reasons:  the  ir-function  strictly  predicts 

improved  performance  only  for  sinusoidal  waveforms.  Square  waves,  as 
discussed  in  the  previous  chapter,  only  represent  a minor  deviation  from 
this  waveform.  In  the  case  of  impulse  forcing  though,  to  still 
administer  the  same  amount  of  drug  over  the  period  in  an  instantaneous 
form  requires  amplitudes  which  are  orders  of  magnitude  higher  than  those 
necessary  for  sinusoidal  forcing.  For  this  reason  another  limitation  of 
the  ir-criterion  becomes  important  here,  the  restriction  that  this  method 
is  valid  only  for  small  amplitudes.  Thus,  the  ir-criterion  probably 
being  of  very  limited  use  in  this  case,  the  only  way  so  far  to  describe 
the  reaction  of  a nonlinear  system  to  impulse  forcing  was  numerical 
integration,  point  by  point,  something  which  is  time  consuming  and 
tedious . 

The  problem  stated  above  was  a major  motivation  for  the  development 
of  the  new  method  presented  in  Section  II. 2.  "Periodic  Impulse  Forcing 
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of  Nonlinear  Systems."  In  the  following  this  method  will  be  applied  to 
cancer  chemotherapy  to  investigate  the  effect  of  drug  injections  on 
normal  and  malignant  cell  population,  to  study  the  effect  of  repeated 
injections,  to  determine  whether  steady  state  or  periodic  operation  is 
superior,  and,  whenever  applicable,  to  determine  the  optimal  forcing 
period . 

The  Effect  of  Drug  Injections  on  Normal  and  Malignant  Tissue 

Assume  a treatment  strategy  where  a drug  is  injected  periodically 
into  the  bloodstream.  This  can  be  described  as  a series  of  impulses 
which  results  in  an  instantaneous  addition  of  a disturbance  m to  the 
state  vector  y_  (an  instantaneous  addition  of  the  dosage  u to  the  plasma 
drug  concentration) . The  modelling  principles  remain  the  same  as  in  the 
previous  chapter.  Together  with  the  specific  cancerous  tissue  to  be 
treated,  the  "limiting  tissue" — the  most  critical  one  of  the  normal  cell 
population — will  be  modelled.  Gompertz  growth  is  assumed  for  both  cell 
types  (equation  4-2)  as  well  as  saturation  effect  towards  the  drug 
(equation  4-3).  Using  first  order  pharmacokinetics  (equation  4-4)  to 
relate  the  drug  dosage  to  the  plasma  drug  concentration  and  treating  the 
body  as  one  compartment,  equations  (4-5)  to  (4-7)  adequately  describe 
this  system.  In  dimensionless  variables  (Appendix  G): 


Y c 
n 


normal  cell  population 


(4-5) 


dt 


Y c 
c 


(4-6) 


dt 


yc  - TVU 


cancer  cell  population 


drug  concentration 


(4-7) 
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To  this  three-dimensional  system  of  equations  the  new  method,  as 
introduced  in  Section  II. 2.,  will  be  applied.  Notice  that  the  equations 
describing  the  cell  populations  are  independent  of  each  other  and  both 
are  dependent  only  on  the  equation  describing  the  drug  instability. 
Therefore  a new,  simpler  two-dimensional  system  can  be  defined 


dy.  Y.c 

dT  = " Vi  ' i = c or  n 


(4-19) 


dc 

dt 


5c 


(4-20) 


Equation  (4-19)  holds  true  for  both  the  normal  cell  population 
(i  = n and  a.  = aQ)  and  for  the  malignant  population  (i  = c and 
ct . = 1).  Having  reduced  the  three-dimensional  system  of  equations  to  a 
two-dimensional  one,  define  now  the  reduced  third  order  Carleman  vector 


This  gives  rise  to  the  third-order  Carleman  matrix  (i  = c,n) 
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(4-22) 


The  vector  of  state  variables  at  the  beginning  of  the  period  is  given  by 
(equation  2-46)  which  in  this  case  becomes 


where 


i-1  * 

Z A.  1 m 
i=0  1 


(4-23) 
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m 


(-m) 


(~m) 

(-m) 


[2] 

[3] 


(4-24) 


with 


(4-25) 


To  administer  the  same  amount  of  drug  over  a period  as  in  the  case  of 
continuous  infusion  (steady-state  administration)  u has  to  be  a function 
of  the  forcing  period  x. 


u = dosage  rate  • x 


(4-26) 
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Now  define  the  A-'s.  According  to  equation  ( 2—  -48) 


A 


o 


(4-27) 


The  reduced  forms  (after  eliminating  redundancies  in  the  monomials  like 


X1X2  and  X2X-j)  of  the  matrices  A^  and  A2  are  as  follows: 
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0 


0 

0 

0 

0 

0 

0 

0 

0 

•3u 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0/ 


(4-28) 


A-  = 


I 0 
0 
0 
0 
0 

0 

0 

3u2 

0 


0 

0 

0 

0 

0 

0 

0 

0 

3u2 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 


(4-29) 


Table  4-1  shows  the  matrix  to  be  inverted  (see  equation  (4-23)).  After 
inversion  it  is  multiplied  by  the  vector  m^  to  find  w,  the  vector  of 
state  variable  at  the  beginning  of  a period. 
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Assume  now  that  the  system  is  subjected  to  a series  of  impulses  and 
that  all  transients  have  subsided.  According  to  equation  (4-2) 


(4-30) 


To  observe  how  the  two  cell  populations  react  to  an  injection  of  drug, 
equation  (4-30)  can  be  solved 


C3T 

vse  So 


i = c ,n 


(4-31 ) 


where 


S = 


(4-32) 


Figure  4-13  shows  the  reaction  of  the  critical  normal  and  the  cancerous 
cell  population  to  an  injection  of  drug  of  the  optimal  steady-state 
dosage  rate  (t  = 1 ) . The  order  of  the  Carleman  approximation 

used  (X,  = 3)  was  found  to  be  very  satisfactory  since  the  results 
obtained  with  this  new  method  are  indistinguishable  (for  most  sets  of 
parameters  within  four  significant  figures)  from  those  obtained  by 
numerical  integration. 


RELATIVE  CELL  POPULATIONS 
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Figure  (4-13) . The  response  of  normal  and  malignant  cells  to 
a single  injection  of  drug  at  the  steady-state 
dosage  ( T=  1.) 
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Table  4-1.  Listing  of  All  Non- Zero  Matrix  Elements  of 


C T 

[e  - I ~ A1  -A2] 

Matrix  element 

Val  ue 

(1,1) 

-a.  t 

e 1 - 1 

0,2) 

x “Qt.  x 

a(e“5x  - e 1 ) 

0,5) 

h,  -25t  "V. 

b(e  - e ) 

0,9) 

ox  -a.x 

o(e"35T  - e 1 ) 

(2,2) 

“6  T 

e - 1 

(3,3) 

-2a.  t 

e 1 - 1 

(3.4) 

(-a. +5 )x  -2a. x 

2a  (e  1 - e 1 ) 

(3,5) 

o "Ot.T 

a2(e  1 - e"5T)2 

(3,8) 

- (a. +25 )x  -2a. x 

2b (e  1 - e 1 ) 

(3,9) 

"a;T  -xt  "a-T  -26x 

2ab(e  - e 6t)(  1 - e ) 

(4,1) 

2u 

(4,4) 

- (a. +5  )x 
e 1 - 1 

(4,5) 

-26x  “(V6)\ 

a(e  - e ) 

(4,9) 

- (a. +6  )x 

b(e  - e ) 

(5,2) 

2u 
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Table  4-1  continued 

Matrix  element 

Val  ue 

(5,5) 

e“26T  - 1 

(6,6) 

~3a.  t 

e 1 - 1 

(6,7) 

-(2a.+5)x  ~3a. x 

3a  (e  1 - e 1 ) 

(6,8) 

-a.x  -a.i  . „ 

3e  x [a(e  1 - e 6t)]2 

(6,9) 

-a.i  . _ 

[-a (e  * - e‘5T)]3 

(7,3) 

3u 

(7,7) 

- (2a. +5  )t 

i 

e 

(7,8) 

-(a.+25)t  ~(2a.+5)x 

2a  (e  1 - e i ) 

(7,9) 

e-5x[a(e"aiT  - e'^)]2 

(8,1) 

-3u2 

(8,4) 

3u 

(8,8) 

- (a. +26  )t 
e x - 1 

(8,9) 

, -(a.+26)t 

ate"35’  - e 1 ) 

(9,2) 

-3u2 

(9,5) 

3u 

(9,9) 

e'36T  - 1 

1 


where 


a 


Y. 

i 


- 6 


b 


c 


i = c ,n 


Optimization  of  Cancer  Chemotherapy 

Again,  the  choice  of  the  proper  performance  measure  is  of  utmost 
importance.  The  ratio  of  the  number  of  cells  at  time  t of  treatment,  as 
represented  by  y^Ct),  divided  by  the  maximum  attainable  number  of  cells 
(without  treatment),  the  performance  measure  is  given  by 


J 


1 [ ,n 

7 > (s 


qe  c)dt 


(4-33) 


q being  the  weight  factor  introduced  in  the  previous  chapter.  This 
performance  measure  has  to  be  expanded  into  a Taylor  series 


J' 


dt 


T 

r w dt 


(4-34) 


where 
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According  to  equation  (2-45) 


J 


1 * r c3  3 


T 

1 + b w 

- -o 


(4-36) 


Table  4-2 


Listing  of  the 

-1  C3T 

c3  [e  3 - I] 


elements  of 


Vector 

Element  Value 


(1) 

(2) 

(3) 

(4) 

(5) 


-a.  t 


1 **  * 
- (e  1 


a. 

l 


D 


(e-<5T 

a. 

l 


-a . t Y. 


e 1 ) + — ~ (e“6T  - 1) 


. -2a. t 

" 4^7  (e  1 - D 

i 


a.  <5 

i 


“(a  +5)t  -2a. t 

■ - — (e  1 
2a. 
i 


- e 1 ) + 


Y. 

i 


2a. (a.  + 6) 
x i 


“(a. +5 )t 
(e  1 


- 1) 


k ~a. t 2 -a . t 

(e  - * 1 ) - hr (e  1 - s'iT>2  * 


4a. 

l 


Y;a  .?rT  “(a.+<5)t 

1 fe  26t  - p 1 ^ 

+ a)  e ) 


2a. (a.  + 5) 

i i 


Y.  Y. 

— (1  + — x-  .J(e~25T 


2a.  5 

l 


2(a.  + 6) 


1) 
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Table  4-2  continued 


Vector 

Element  Value 


(6) 

(7) 

(8) 


1 "3V 

TfoT  (e  1} 

i 

-(2a.+5)x  -3a. x Y. 

a , i i . i 

6a.  6 6 6a. (2a.  + 5) 

1 XL 


- (2a. +5 )t 


(e 


- 1) 


■a.i  -a.i 


b , i 

2^7  (e 

l 


-<c,  ♦ 2«)t  -2Vi  3e  -i-fa(e  -i'  . e-it;12 


- e ) - 

aY.d  -(a.+26)x  -(2a.+5)x 


1 8a 


3a ; 


i / i 

— (e 


- e 


) 


Y. 

i 


(a.  + 26 )a.  2 3 

i i 


. Y.d  - (a. +26  )t 

, 1 i , , i 

(~  + -—)  (e 


1) 


(9) 


„ -a.t  , -a.T  . -a.T  or 

c , ~36x  i ..  ab  , x -5tw  i -26tn 

(e  - e ) - -r — (e  - e )(e  - e ) + 

a.  2a. 

x x 


Y.b 

i 


2a. (a.  + 6 ) 
i x 


-(a.+6)x 

(e'3iT  - a 1 ) 


1 


fa(e  1 - e'6T)13  

18a.  6a. (2a.  + 5 ) 

l ii 


-a.  x 


[a (e  1 - e~6T)]2 


Y.  a 

i 


a* 


Y. 

1 


(a.  + 26)a.  2 3(2a.  + 6) 

ii  i 


_ -(a.+26)t 

/ ~36t  i . 

(e  - e ) 


Y.  Y.  . 

— i_  [i  + ( 1 

3a, 6 L 2 a.  + 


-)  + 


Y2 

1 


26  a + 5 3(2a.  + 5)(a.  + 26) 

l ii 


<e-3ST  - 1) 


a,  b,  c are  as  defined  in  Table  4-1. 


109 


Figure  4-14  shows  that  for  Yn  = 350,  Yq  = 401 , 6 = 6.67,  a = 5,  q = .2 

periodic  treatment  indeed  does  yield  better  results  than  continuous 

treatment  as  expected  from  the  ir-criterion.  Furthermore  it  is  obvious 

that  the  new  method  used  in  this  chapter  is  in  excellent  agreement  with 

the  mathematical  simulations  (integrations)  in  predicting  the  optimal 

length  of  the  treatment  cycle  (t  . = 2.2)  while  the  predictions  of 

opt 

the  iT-criterion  are  le3S  accurate  in  this  case  (t  . = 1.6). 

opt 

For  a different  set  of  parameters  however,  (Y^  = 240,  Yc  = 401,  6 = 

66.7,  a = 8,  q = .5)  the  new  method  fails  to  predict  the  optimal 

treatment  frequency.  It  still  yields  the  vector  w^,  the  cell 

populations  at  the  beginning  of  each  treatment  cycle  with  highest 

accuracy,  but  in  this  case  the  treatment  is  so  successful  that  |yc|  > 

1.  This  of  course  means  that  the  Taylor  series  expansion  used  in  the 

performance  measure  is  no  longer  valid  and  thus  the  optimal  length  of 

the  treatment  cycle  cannot  be  predicted.  However  in  this  case 

the  ir-criterion  succeeds  in  predicting  the  optimal  period  (t  t = 0.87 

opt 

versus  0.85  obtained  by  integration). 

The  new  method  has  particular  value  for  all  those  cases,  where 
the  ir-criterion  is  inapplicable  because  it  does  not  require  (like 
the  ir-criterion  does)  that  the  performance  measure  J exhibits  a maximum 
with  respect  to  the  control  variable  u.  If  instead  of  the  performance 
measure  (4-33)  one  chooses  to  maximize  the  ratio  of  the  relative  cell 
populations,  as  might  be  indicated  for  some  types  of  cancer,  a linear 
performance  measure  arises. 


x /x  ,<*> 


J = - f (in  — 7 )dt  = - J (y.-y-Mt 

T J V X /X  T J 12 

o c c o 


(4-37) 
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Figure  (4-14).  The  prediction  of  the  optimal  cycling  period 
by:  a)  numerical  integration 

b)  the  new  method 

c)  the  ir-criterion 


Ill 


No  optimal  steady-state  value  exists  for  the-  control  variable  and 
therefore  the  ir-criterion,  which  compares  periodic  operation  to  the 
optimal  steady-state  operation,  cannot  be  applied.  For  the  performance 
measure  as  defined  by  equation  (4-38) 


1 fT  T r^  -1  ^3t  t 

- J r w dt  = — [e  J - I]  w_  = b1  w. 


-o  - -o 


(4-38) 


where  r is  given  by 


r = 


(4-39) 


In  this  case  the  vector  _b  in  equation  (4-38)  is  given  by 


/ 1 "V 

- 1 (e  1 - 1) 


®i  . -a.T  Y. 

- J-  (e'5T  - e 1 ) . -i  (e-'T  - 1) 
ot.  a.  <5 

i l 

0 


K rt.  -a.x  Y. 

. £.  (e"25'  - e 1 ) - ^ (e'2ST  - 1) 
a.  2a. 6 

i i 

0 
0 
0 

1 1 


( 4-40) 
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Figure  4-15  3hows  the  dependence  of  the  linear  performance  measure  J on 

the  cycling  period  t for  a certain  set  of  parameters  (y  = 240,  Y = 

n c 

401,  6 = 66.7,  a = 8,  q = 5).  Three  dosages  were  investigated,  and  the 
results  show  clearly,  that  continuous  infusion  (steady-state  treatment) 
is  superior  to  any  kind  of  periodic  regimen.  The  higher  the  dosage,  the 
higher  is  the  performance  measure.  The  optimal  dosage  rate  therefore  is 
the  highest  one  allowed  by  toxicity  criteria.  In  most  cases  continuous 
treatment  is  rather  inconvenient,  because  it  may  force  the  patient  to 
stay  in  the  hospital  for  an  extended  amount  of  time.  Alternatively,  if 
this  form  of  schedule  is  not  possible,  the  treatment  cycles  should  be 
kept  as  short  and  as  frequent  as  possible  and  the  dosage  the  highest 
possible  one. 

Conclusions 

In  the  case  of  impulse  forcing,  because  of  the  strong  deviation 
from  the  sinusoidal  wave  form  and  the  high  amplitudes  necessary, 
the  u-criterion  does  not  always  succeed  in  predicting  the  optimal  length 
of  the  treatment  cycle.  The  new  method,  on  the  other  hand,  is  capable 
of  predicting  the  optimal  treatment  strategy  accurately,  as  long  as  the 
state  variables  are  defined  in  such  a way  that  they  remain  in  the  range 
where  the  Taylor  Series  expansion  is  valid. 

For  a performance  measure  which  maximizes  the  difference  in  the 
cell  populations,  weighed  by  a weighing  factor,  a periodic  schedule  is 
superior  to  continuous  treatment  and  an  optimal  cycle-length  can  be 
found.  When  a linear  performance  measure  is  employed,  continuous 
treatment  is  preferable.  If  this  is  not  possible,  the  treatment  cycles 
should  be  kept  as  short  and  the  dosage  rates  as  high  as  possible. 


PERFORMANCE  MEASURE 
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Figure  (4-15).  The  effect  of  the  dosage  on  the  (linear) 
performance  measure 


CHAPTER  V 

SUMMARY  AND  CONCLUSIONS 


Over  the  past  few  years  it  has  been  realized  that  periodic 
operation  of  biochemical  reactors  can  be  superior  to  steady-state 
operation.  Although  some  workers  had  tried  to  assess  the  effect 
experimentally  in  the  past,  no  systematical  investigation  of  the  effect 
of  periodic  operation  had  ever  been  undertaken. 

In  a very  different  scientific  field,  cancer  chemotherapy,  it  had 
been  disputed  for  many  years  whether  cyclic  or  continuous  treatment  of 
cancer  will  yield  better  results.  In  this  dissertation  the  attempt  has 
been  made  to  answer  some  of  these  questions. 

There  are  two  possible  approaches  to  determining  the  optimal 
periodic  conditions:  an  ad  hoc  trial-and-error  procedure  or  a 

systematical  search  utilizing  mathematical  modeling  as  a guidance 
tool.  Based  on  mathematical  process  models,  methods  of  optimal  control 
theory  can  be  applied  to  answer  the  question  of  what  the  optimal  period, 
amplitude  and  wave  form  are. 

Among  the  already  established  methods  perhaps  the  most  well  known 
is  the  ir-criterion.  It  has  proven  useful  in  many  cases  for  finding  the 
optimal  periodic  operating  conditions,  but  does  not  perform  quite  as 
well  in  cases  where  the  wave  form  deviates  strongly  from  a sinusoid  or 
where  amplitudes  are  too  high. 


114 


115 


For  this  reason  a new  method  "periodic  impulse  forcing  of  nonlinear 
systems"  has  been  introduced  in  this  dissertation  specifically  for  the 
study  of  periodic  impulse  forcing. 

The  ir-criterion  was  applied  to  predict  optimal  periodic  conditions 
for  a chemostat  undergoing  square  wave  forcing  in  the  dilution  rate.  It 
was  concluded  that  enhanced  productivity  of  biomass  can  restilt  for 
certain  cycling  frequencies.  This  is  a consequence  of  the  fact  that 
there  is  a time-lag  in  the  culture  adaptation  to  the  new  growth  rate 
when  a disturbance  in  introduced.  The  predictions  obtained  with  methods 
of  optimal  control  theory  can  serve  as  a guide  for  finding  the  actual 
experimental  operating  conditions,  which  maximize  the  performance 
measure.  In  order  to  be  able  to  describe  the  experimental  results 
qualitatively  it  was  indeed  necessary  to  include  a time-lag  in  the 
unstructured  model  describing  the  response  to  a step  change  in  the 
dilution  rate.  It  was  found  that  for  a step  down  the  culture  continues 
to  grow  according  to  the  initial  growth  rate  for  a certain  time,  while 
when  the  dilution  rate  was  stepped  up,  growth  even  stalled  during  this 
lag-time . 

The  same  effects  were  observed  when  the  system  was  subjected  to 
square  wave  variations  in  the  dilution  rate.  A considerable  improvement 
in  the  reactor  performance  was  obtained  for  a cycling  frequency  of  three 
hours.  A possible  explanation  for  this  phenomenon  is  that  in  this  case, 
a high  (exponential)  growth  rate  coincides  with  a high  concentration  of 
biomass  at  a time  when  there  is  ample  substrate  present. 

Just  as  the  growth  of  cells  in  a reactor  can  be  maximized  by  a 
certain  "feeding  schedule,"  the  growth  of  unwanted  cells,  i.e.  cancer 
cells,  can  be  minimized  by  a certain  treatment  schedule.  The  limiting 
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factor  in  cancer  chemotherapy  is  the  toxicity  of  the  drug  to  the  normal, 
homeostatically  controlled,  cell  population.  Thus,  an  optimization 
problem  exists  between  killing  the  maximum  number  of  cancer  cells  while 
doing  as  little  damage  as  possible  to  the  benign  cell  population. 

The  ir-criterion  helped  in  finding  an  answer  to  the  following 
questions:  (1)  Is  the  drug  appropriate?  (2)  Will  continuous  or  cyclic 

administration  give  better  results?  (3)  How  much  shall  be  administered 
and  in  which  time  intervals?  Although  the  mathematical  model  used  was 
rather  rudimentary,  it  was  based  on  actual  patient  data  and  the 
resulting  optimal  treatment  schedule  is  among  those  used  in  clinical 
practice. 

If  the  point  of  interest  is  the  reaction  of  normal  and  malignant 
cells  to  a single  injection  of  drug  or  a series  of  injections,  the 
predictions  of  the  ir-criterion  were  found  to  be  less  accurate  because 
some  of  its  limitations  became  relevant  in  this  case.  This  was  a major 
motivation  for  the  development  of  the  new  method  "periodic  impulse 
forcing  of  nonlinear  systems."  This  method  has  been  found  to  be  very 
successful  in  predicting  the  optimal  periodic  conditions,  as  long  as  the 
state  variables  are  defined  in  such  a way  that  they  remian  in  the  range, 
where  the  Taylor  Series  expansion  is  valid. 

The  quality  of  any  prediction  made  by  methods  of  optimal  control 
theory  depends  to  a great  extent  on  the  quality  of  the  mathematical 
model  used  to  describe  the  system  in  question.  This  in  turn  depends  on 
the  quality  of  the  experimental/patient  data.  The  more  complete  the 
available  data,  the  better  is  the  mathematical  model  and  the  better  are 
the  predictions  about  the  optimal  periodic  operating  conditions. 
Therefore,  in  optimizing  any  process  it  frequently  is  worth  the  effort 
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to  conduct  experiments  designed  especially  to  obtain  the  best  possible 
dynamic  model  in  an  effort  to  obtain  predictions,  which  are  not  only 
qualitatively  true,  but  also  quantitatively  good. 


APPENDIX  A 

OPTIMAL  PERIODIC  SQUARE  WAVE  FORCING: 

A NEW  METHOD  (LYBERATOS  AND  SVORONOS  [24]) 

Consider  a process  described  by 

X = F(x,u);  F(0,0)  = 0 (A-1) 

where  £(x_,  u)  is  analytic  in  £ at  £ for  all  values  of  u in  the 
admissible  range.  £ and  u are  the  3tate  vector  and  scalar  input 
deviations  from  their  optimal  steady-state,  respectively.  The  objective 
is  to  find  the  periodic  pulse  input  described  by 


«-eT-fc 


o te  [nT,  (n+e )T] 

c5 

P = — f te  [(n+e)T,  (n+1  )T] 


(A-2) 
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that  maximizes  the  time-averaged  performance  measure 


J = ^ / P(x,  u)dt  (A-3) 

0 

which,  in  Carleman  coordinates,  takes  the  form 

J = ? tr0(u)  + r:T(u)w(t)]dt  ( a— 4 ) 

w(t)  being  the  state  vector  (in  Carleman  coordinates)  as  introduced  in 
Chapter  II. 


Theorem:  The  performance  measure  is  given  by 


J(T,  e,  5)  =[rT(p ) S(p)"1  - rT(6)  1- (I~R) iLMl  ^'D) 

T 


[S(p)  ^(p)  - S ( 6 )_1z((5 ) ] 


-1 


(A-5) 


S(p)  ^(p)  - rT(6)  S( 6 )" 1 z (5  ) ] 


T -1 

- r (p)  S(p)  z(p)  + rQ(6)e  + r (p)(l-e) 

where  R = exp  [s(p)(1-e)T]  and 
D = exp  [s(6)eT] 

and  3(*)  is  the  Carleman  matrix. 


APPENDIX  B 

DERIVATION  OF  THE  MATRICES  Ai 
(CHAPTER  II. 2.) 


This  appendix  will  illustrate  more  explicitly  how  equation  (2-49) 
has  to  be  used.  In  the  following  the  matrices  AQ  - Aj_1  will  be  derived 
for  an  n-dimensional  system  and  the  4th  order  of  the  Carleman 
approximation  (£=4) 

by  definition:  A«  = I 


o 


n x n 


% 


A.  = 
i 


0 


l 

n x n 


therefore 
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A,  = 


0 

-n  x n 


E I®  (-m) 


n3  x n 


2 4 

n x n 


0 


n x n 

-°2  2 
n x n 

E IL^J®  (-m) 


- 4 2 

n x n 


0 


n x rr 


- 2 3 

n x n 


n3  x n3 


n x n 


■ 2 4 

n x n 


-3  4 

n x n 


* m 

EIUJ®(-m)  o ^ 

n x 


0 


A2  = 


n x n 


32 

n x n 


E*I  ® (-m)[2] 


2 4 

n x n 


n x n 


-2  2 
n x n 


n3  x n2 


0 


n x n' 


-°  2 3 
n x n 

~n3  x n3 


Z*IC2]®  (-m)C2]  0 


4 3 

n x n 


0 


n x n 


-2  4 

n x n 


- 3 4 

n x n 


- 4 4 

n x n 


#3  * 


0 


•n  x n 


- 2 

n x n 


23 

n x n 


0 


n x n 


- 2 2 
n x n 


V * n2 


El®  (-m)^  0 


-4  2 

n x n 


n x n' 


■ 2 3 

n x n 


n3  x n3 


-4  3 

n x n 


0 


n x n 


-2  4 

n x n 


; 3 4 

n x n 


■4 

n x n 


Ai  always  is  a square  matrix  with  the  dimensions  of 


l nk. 
k=1 


APPENDIX  C 

DERIVATION  OF  RESULTS  IN  TABLE  (2-1) 


This  appendix  will  illustrate  how  the  results  in  Table  2-1  were 
obtained.  They  will  be  derived  for  the  third-order  Carleman 
approximation.  In  dimensionless  variables,  the  system  described  by  eq. 
(2-56)  around  x=1  is 


x = - x - 


2 

x 


(C-1) 


The  third-order  Carleman  vector  is  defined  as 


w 

-o 


x. 

x , 

1 

1 

Y 

XC2] 

X2 

X1 

Y 

x[3] 

3 

X1 

(C-2) 


2x2  - 2x 
- 3x 


(only  terms  up  to  3r“d  order 
were  retained) 


Therefore, 


-1  -1  0 


-2 


-2 

-3 
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0 

0 


0 


(C-3) 
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From  eq.  (2-41) 


w = 
-o 


* 

Z A.  ] m 
i=0  1 


(C-4) 


where 


0 0 0 
-2m  0 0 

0 -3m  0 


and 


(C-5) 


(C-6) 


so 


C3t 


- 1 - Ar  *2l 


-1  * 
m 


(C-7) 


-T 


e 2x  - e~T  e~3x  - 2e~2x  + e”x 


■2t 


2e"3x  - 2e"2x 


“3t 


(C-8) 


Therefore  the  matrix  to  invert  is 


E = 


e_x  - 1 e"2x  - e“x 


2m 


-3m 


e_2x  - 1 


e'3x  - 2e'2x  + e'x 


2e“3x  - 2e"2x 


-3t 


(C—  9 ) 


3m 


1 
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For  small  values  of  x the  matrix  E may  become  nearly  singular,  which  can 
lead  to  numerical  errors  during  inversion.  In  this  case  matrix  E needs 

to  be  evaluated  and  inverted  numerically,  and  numerically  multiplied  by 

£ 

the  vector  m . 


(C- 10) 


where 


(-m) 

(-m) 

(-m) 


[2] 

[3] 


(C-11) 


This  can  be  done  using  a maximum  pivot  strategy.  The  performance 
measure  is  given  by 


(C- 12) 


where 


1 

0 

0 


(C- 13) 


Combining  eq.  (C-10)  and  (C-11)  or  (C— 12)  an  (C-13).  respectively,  leads 
to  the  results  in  Table  2-1. 


APPENDIX  D 

*(“)  FOR  A MONOD  MODEL  WITH  DELAY  AS  GIVEN  BY 
EQUATIONS  (3-11)  - (3-13) 


/ \ 2 r 4 2 2 2 “*9  ^ ^ o o'*  o 

ir(o))  = - ia)  F GJ  + u FG (D  rGK  + DfVk  + D2FG-  DGa2  + 


DFGK  a 
s 


where 


- FGK  a2)  - D2FG2a2(D  + FK  ) } 

s 3 i 


6 2*2  4 9^9  ^ ~ 9 

M = w F Ks+  0)  F K [(a+D+FG)  - 2a(D+FG)  ] 


2 2a2a  a " 2 
+ a)  F K a[a(D+FG ) - 2DFG (a+D+FG ) ] 


~2  4 2 "2  "2 
+ D F G K a 

3 


F = 


1-D 


and 


G = 1-D-DK 


In  the  above  D is  the  optimum  steady-state  dilution  rate, 
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APPENDIX  E 

tt(u>)  FOR  MONOD  MODEL  WITH  DELAY  AS  GIVEN  BY 
EQUATIONS  (3-11),  (3-12)  AND  (3-21) 


it(u>) 


4D2  - 4D2eos(u)T)y  2K2x^ 
m s 

2 4 
ADY  F 


2Ey  K x w oosojt 

m 3 

ADYF2 

[2D-4DOOS  (u>t  ) ]DymKsX2  2DKgymX2wsin  Ut  ) 
AYF2  AYF2 


2 2 2 
2Dx  (D  + on  ) 

AYF 


where  F = K + s , 

3 


3 


K D 


D 


X = Y (S°  - 3) 


and 


A = ( 


DK  y xcos(u)t) 

3 m 

YF2 


2 2 DK  y xsin(ux)  0 

u2)2  + ( + Du  _ _ ?.  m )2 


YF 


YF 


In  the  above  D is  the  optimal  steady-state  dilation  rate. 
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APPENDIX  F 

ir(u)  FOR  THE  WILLIAMS  MODEL 


= jjf  [-32w  a2(1  +a)^  + 8<D2(4a5  + 160^+  24a^+  18a2+  7a  + 1) 


+ 8a^  + 1 6a2+  10a  + 2] 


where 


M = (-4w2a2+  1)  • [[(1  + 2a)  - 4(1  + a)2aj2]2+  16(1  -*-a)4uj2] 
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APPENDIX  G 

DERIVATION  OF  THE  DIMENSIONLESS  MODEL  FOR 
CANCER  CHEMOTHERAPY  (EQUATIONS  (4-5)  TO  ((4-8)) 


Given  the  Gompentz  equation: 


dUnx. 

1 

dt 


A 


0 


x. 

- a.  In  1 - 
i 

x. 

1 ,o 


(G-1) 


defining 


A 


o 


a.  £n 
i 


x. 

1 ,o 


(G-2) 


Substituting  ( A2)  into  (A1) 


dx.  x. 

- x a S,n  — (G-3) 

dt  1 x. 

i 


defining  dimensionless  variables  as  follows: 
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y = fcn  — 
n x 


n ,® 


k 

Y„  = - 
n a 

o 


e 


c 

k 


t 


= a t 
o 


and  substituting  into  (G-4  - G-6) 


dx 

n 

/\ 

dt 


a x £,n 
n n 


k c 
n 


k + c 


(G—  4 ) 


dx 

e 

dt 


a x £n 
c c 


k c 
o 

* x 

. 0 
k + c 


(G-5) 


d~ 

-^  = u - 6c  (G-6) 

dt 


gives  the  dimensionless  model 


dt 


ayn  - 


Y C 
G 

1 + G 


(G-7) 


dyc  Yo° 

dt  - yc  1 + o (G-8) 


do 

dt 


u 


5c 


APPENDIX  H 

CALCULATION  OF  THE  OPTIMAL  TREATMENT  PROTOCOL 
FOR  CANCER  CHEMOTHERAPY  (CHAPTER  IV) 

Keeping  the  total  amount  of  drug  constant  oven  the  length  of  a 
treatment  cycle,  e is  defined  as  the  fraction  of  the  cycle  the  patient 
receives  a drug  (Fig.  4-11). 

The  dosage  w given  during  this  period  is  u /e  (so  that  equation 

(4-15)  is  satisfied).  For  small  values  of  e,  u might  become  higher  than 

the  highest  physiological  allowable  dosage  Then  u = umax 

and  e . = 

mm  u 

max 

An  example  shall  illustrate  this  method:  For  the  following 

parameter  choices: 


a = 15 

Y 

= 401 

c 

6 = 66.7 

Y 

= 160 

n 


the  iT-criterion  yields  an  optimal  treatment  cycle  of  35  days.  Numerical 
integrations  show  that  e ^ = 0.08  which  translates  into  3 days  of 
treatment  and  32  days  without  medication. 
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